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ABSTRACT 

The  classical  Levinson-Durbin  linear  prediction  formulas  for  real  valued  input  sequences 
are  examined  and  compared  to  the  recently  proposed  split-Levinson  formulas.  Both  the 
autoregressive  linear  predictor  model  and  the  adaptive  lattice  model  are  used  to  formu- 
late the  new  split-Levinson  algorithms.  A  brief  introduction  to  the  theory  of  symmetric 
polynomials  is  presented  to  form  the  basis  of  the  new  algorithms.  Computer  simulations 
are  used  to  test  and  compare  the  computational  accuracy  of  the  new  algorithms  for  AR 
filter  coefficient  estimation,  parameter  estimation  for  a  moving  average  process,  and 
spectral  estimation  of  sinusoids  in  white  noise.  Research  results  indicate  that  the  new 
algorithms  reduce  the  number  of  real  multiplications  required  for  a  k"'  order  AR  filter 
problem  by  one-half,  and  they  are  applicable  to  both  the  extended  Prony  method  of 
spectral  estimation  and  the  estimation  of  moving  average  parameters. 
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I.     INTRODUCTION 

A.      OBJECTIVE 

The  classical  Levinson  algorithm  is  known  to  provide  solutions  to  real  valued,  linear 
systems  involving  Toeplitz  structures.  The  computational  cost  for  these  solutions,  is 
known  to  be  0{k:).  where  k  indicates  the  filter  order.  It  has  recently  been  proposed  that 
the  classical  algorithm  may  be  transformed  into  2  simpler  algorithms,  using  the  theory 
of  symmetric  polynomials,  and  that  either  of  these  algorithms  can  be  used  to  to  solve  for 
the  predictor  polynomial  of  order  k  at  a  reduced  computational  cost.  [Ref.  1:  p.  470] 

These  new  algorithms  are  termed  the  split-Levinson  algorithms  because  their  basis 
is  formed  from  the  concept  of  symmetric  polynomials.  These  are  not  new  in  theory,  but 
the  application  of  the  process  to  linear  prediction  is  a  new  concept.  Symmetric 
polynomials  are  based  on  the  Barlett  Bisection  Theorem  [Ref.  2  :  pp.  1074-1076],  where 
a  system  that  possesses  symmetry  about  a  point,  such  as  a  Toeplitz  matrix,  can  be  de- 
composed into  a  symmetric  and  an  antisymmetric  part.  The  unique  point  of  the  theory 
is  that  either  part  may  be  used  to  solve  the  problem,  or  a  combination  of  both  parts  can 
also  be  used  in  the  solution.  During  our  research  we  shall  only  consider  real  data  se- 
quences. 

The  split-Levinson  case  also  has  a  lattice  structure  as  the  classical  case.  However, 
as  will  be  shown,  the  structure  of  this  lattice  shows  little  resemblance  to  its  classical 
counterpart.  A  derivation  of  a  revised  split  lattice  structure,  and  its  recursive  algorithm 
was  attempted  in  order  to  represent  the  split  lattice  in  a  form  similar  to  that  of  the 
classical  structure.  Although  unsuccessful,  the  derivation  procedure  is  presented  for 
subject  matter  continuity. 

Computer  programs  have  been  written  to  implement  the  new  algorithms  and  com- 
pare them  to  the  classical  algorithms.  Additionally,  computer  programs  are  included  to 
apply  the  split-Levinson  algorithm  for  two  cases,  where  the  computational  efficiency  of 
the  new  algorithm  could  be  of  substantial  benefit.  These  cases  include  the  Moving  Av- 
erage (MA)  problem,  where  the  parameters  of  a  MA  model  must  be  determined  from  the 
given  data,  and  an  extension  of  the  Prony  method  of  spectral  estimation,  where  a  least 
squares  estimation  of  the  presence  of  sinusoids  in  white  noise  is  made  from  the  output 
data  sequence. 


This  thesis  compares  the  classical  and  lattice  structures  of  the  Levinson  recursion 
formula  given  in  [Ref.  3:  pp.  145-167],  and  examines  not  only  the  formulation  of  the 
recursion  formulas  for  these  algorithms,  but  also  the  complexity  o^  the  computations 
and  the  resulting  structure  of  each  of  the  algorithms. 

B.     THESIS  ORGANIZATION 

The  structure  of  the  thesis  is  divided  into  4  chapters,  including  the  Introduction.  In 
Chapter  II  we  will  review  the  classical  Levinson  algorithms.  In  the  first  case,  the  algo- 
rithm is  obtained  using  the  autocorrelation  function  of  the  input  sequence,  and  in  the 
second  case,  it  is  obtained  using  the  forward  and  backward  error  vectors  of  the  input 
sequence.  In  each  case  we  shall  establish  the  number  of  real  multiplications  required  to 
complete  a  k-th  order  recursion  of  the  respective  algorithm.  As  stated,  the  ultimate  goal 
is  to  establish  the  computational  efficiency  of  the  split-Levinson  algorithm  over  the 
classical  Levinson  algorithm.  Chapter  III  deals  with  the  derivation  of  the  split-Levinson 
algorithms  preceded  by  an  introductory  section  on  symmetric  polynomials.  As  in  Chap- 
ter II,  both  the  autocorrelation  function  and  the  lattice  algorithms  will  be  developed. 
In  addition,  a  comparison  between  the  computational  cost  of  the  Levinson  and  split- 
Levinson  algorithms  and  an  attempt  to  define  the  split  lattice  structure  in  terms  similar 
to  the  Levinson  based  lattice  are  presented. 

In  the  final  chapter  two  practical  applications  of  the  split-Levinson  algorithm  are 
investigated.  These  are:  (1)  the  MA  parameter  estimation  problem,  and  (2)  the  extended 
Prony  method.  In  case  (1),  the  Levinson  recursion  used  to  determine  a  predictor  coeffi- 
cient vector  is  replaced  by  the  split-Levinson  algorithm.  A  comparison  between  the  test 
coefficients  and  the  computed  coefficients  is  presented.  In  case  (2),  an  estimation  of 
sinusoids  in  white  noise  is  performed.  Additionally,  overall  conclusions  of  the  research 
as  well  as  proposed  topics  for  continued  thesis  research  are  presented. 


II.     THE  CLASSICAL  LEVINSON  ALGORITHMS 


The  importance  of  the  Lcvinson  algorithm  in  linear  prediction  theory  is  well  known. 
The  reason  to  present  the  algorithm  in  us  two  forms  is  twofold:  (1)  to  present  certain 
definitions  that  will  be  required  later  in  the  development  of  the  split- Le Vinson  algo- 
rithms, and  (2)  to  detail  the  computational  complexity  of  the  Levinson  algorithm  for 
comparison  purposes  to  the  split-Levinson  versions  of  the  algorithm.  In  the  context  of 
our  discussion  within  this  thesis,  we  shall  confine  ourselves  to  the  study  of  autoregressive 
modeling  problems  of  real  sequences  as  in  Figure  1.   [Ref.  3:  p.  152] 


Figure   1.        Autoresressive  Model. 


We  know  from  linear  prediction  theory  the  augmented  normal  equations  given  by, 

Rx*  =  rxy  (2.1) 

can  optimally  be  solved  by  the  Levinson  algorithm,  and  that  this  algorithm  can  be  im- 
plemented with  either  the  autocorrelation  function  or  the  forward  and  backward  pre- 
diction error  vectors  of  the  input  sequence  [Ref.  3:  pp.  152-170]. 

A.      THE  LEVINSON  ALGORITHM 

In  order  to  examine  the  computational  complexity  of  the  Levinson  recursion,  it  is 
necessary  to  formulate  the  recursive  algorithm,  and  to  determine  the  number  of  real 
multiplications  and  additions  required  to  complete  the  algorithm.  First,  construct  a 
Toeplitz  matrix  from  the  sequence  s(t),  of  length  N,  defined  as  Rk  =  [/?,_,  :0  <  /,  j<  fc], 
where  the  elements  of  the  matrix  are  the  autocorrelation  lags  given  by  [Ref.  2:  p.  646] 

(V- 1  -i 
Ri=  V-i-i    Z    *W«'+0  (2-2) 

i'=0 

then,  the  predictor  vector  a  can  be  determined  as  a  solution  to  the  system  defined  by  the 
matrix  equations 

:k*:m-i>4,o,o of  (2.3) 

where  ak  is  the  prediction  error  norm,  and  is  defined  as 


°  k 


It  is  recalled  from  linear  prediction  theory,  that  given  a  positive-definite  matrix  Rk 
of  order  k+  1,  the  kth  order  a  coefficient  vector  can  be  computed  recursively  from  the 
nested  Toeplitz  submatrices,  and  their  respective  successive  predictor  vectors,  a.  The 
well-known  Levinson  recursion  formula  is  this  solution,  and  has  the  form 

akl  =  ak_]i  +  pkak_]k_i         /  =  0,1,2,...,/:  (2.5) 

with  the  conditions  that  aK0=  1,  and  ak_Uk  =  0.  The  parameters,  pk  =  akj!,  are  called  re- 
flection coefficients,  also  PARCOR  coefficients,  because  they  represent  the  partial  cor- 


relation  between  the  zero-th  and  the  k-th  element  of  the  prediction  vector  with  the  effect 
of  all  the  intermediate  elements  removed. [Ref.  4:  p.  53] 

To  construct  the  Levinson  recursion  we  must  use  the  prediction  error  norm  re- 
lationship 

*A  =  (1-Pjfc)*fc-1  f2-6) 

and  the  identity 

fe-i 

ok_xpk  =  -  l^  Rk-{ak_xl  (2.7) 

(=0 

to  define  the  recursion  variables.  Consider  the  following  definition  as  it  applies  to  the 
Levinson  recursion  [Ref.  1:  p.  472]. 

k-\ 

lk  =  ~  2_j    Rk-iak-l,i  (2  8) 

(=0  v    '    ' 

=  ak-\pk 

and  solving  for  pk  from  Eq.(2.8),  we  have 

pk-^tr  (2.9) 

The  error  norm  ak  can  be  written  in  terms  of  the  the  normalizing  term  Xk  by  rewriting 
Eq.  (2.6),  and  making  a  substitution  from  Eq.  (2.9) 

CT*  =  (1  ~  Pk)°k-\ 

=  <?/<_!-  Pk(°k-\Pk)  I2-10) 

=  ok-x  -  pkh 

Combining  Eqs.  (2.5),  (2.8),  (2.9),  and  (2.10),  we  have  the  basis  for  the  Levinson  algo- 
rithm, and  it  is  summarized  in  Table  2  of  Appendix  A. 

B.      LEVINSON  LATTICE  REALIZATION 

If  we  are  given  a  real  sequence  of  signal  values  s(0).s(l),  ...,  s(N-l),  and  it  is  known 
that  s(t)  =  0,  for  —  1  >  I  and  /  >  N,  then  for  the  linear  prediction  problem  of  order  k.  we 
find  it  necessary  to  find  a  set  of  real  numbers  ak0.  aku  ... ,  akk  that  will  minimize  the  for- 
ward and  backward  prediction  errorvectors  using  a  linear  combination  of  the  past  signal 


vectors.  If  we  call  the  forward  prediction  error  vector/^;)  and  the  backward  error  vector 
bk{t),  and  define  them  in  terms  of  the  akl  coefficients.  [Ref.  2:  p.  646]  we  have 


Aw=y%  s{t~f)  (2.H) 


(=0 

then  it  turns  out  that  the  same  real  numbers,  ak,  ,  will  provide  the  solution  to  either  of 
the  forward  or  backward  prediction  problems,  (i.e.,  minimize  the  squared  Euclidean 
norm  of  both/A  and  bk). 

Let  ak  be  defined  as  the  squared  norm,  that  is 

a,  =  ||y2  =  f[f,  (2.13) 

From  Eqs.  (2.11)  and  (2.12)  forming  the  first  three  terms  of  each  error  vector  we  have 
the  following, 

MO)  =  akos(0)  +  akls{  -1)  +  ...  +  akks(  -k)  (2.14) 

fk(l)  =  akos(\)  +  akls(0)  +  ak2s(  -1)  +  ....  +  akks{\  -  k)  (2.15) 

fk(2)  =  akos(2)  +  akls{\)  +  ak2s(0)  +  ...  +  akks(2  -  k)  (2.16) 


and, 


bk(0)  =  akks(0)  +  akJi_lS(  -1)  +  akJ(_2s{  -2)  +  ...  +  W(  -k)  (2.17) 

bk{\)  =  akks{\)  +  akk_,s(0)  +  akJc_2s(  -1)  +...  +  akos(i  -  k)  (2.18) 

**(2)  =  auM)  +  "kjc-AV  +  "kji-iW)  +  -  +  W(2  -  k)  (2.19) 


If  we  examine  the  elements  of  these  two  vectors  we  can  see  they  are  related  in  that 
each  can  be  derived  from  the  other  by  reversing  the  order  of  the  a  coefficients.  If  we  form 
the  Euclidean  norm  of  each  vector,  ||fj|  and  ||bj|,  we  see  that  the  k-th  predictor  vector 
\am,  akx akky  minimizes  the  error  norm,  and  \[fk\\  =  \\bk\\. 

From  [Ref.  3:  pp.  156-157],  we  can  use  the  Levinson  algorithm  to  define  the 
recursion  formula  for  the  forward  and  backward  prediction  errors  given  by 

fk(t)=fk-l(t)  +  Pkbk_l(t-  1), 

bk(t)  =  plfk_lU)  +  bk_](t-l)  U-J) 

If  we  let  the  following  definition  apply  to  the  lattice  version  of  the  Levinson  algo- 
rithm 

\+k-2 

h  =  ok_xPk  =  -    Y^  fk_,{t)bk_x{t)  (2.21) 

r=i 

then  using  Eqs.  (2.9),  (2.10),  (2.20),  and  (2.21)  we  can  summarize  the  Levinson  lattice 
algorithm  in  Table  3  of  Appendix  A. 

Even  though  the  lattice  algorithm  is  implemented  directly  from  the  data  samples,  its 
computer  implementation  will  be  more  complex  because  of  the  vector  manipulations 
that  must  occur  in  each  iteration.  The  lattice  structure  defined  by  Eq.  (2.20)  is  shown  in 
Figure  2. 

In  summary,  we  discussed  the  Levinson  algorithm  which  uses  the  autocorrelation 
elements  in  its  recursion  and  the  related  lattice  structure  which  uses  the  input  data  di- 
rectly in  its  formulation.  In  terms  of  the  computational  complexity,  both  algorithms  re- 
quire real  multiplications  of  the  order  k2,  as  detailed  in  Tables  2  and  3  of  Appendix  A, 
in  order  to  realize  a  k'h  order  predictor  filter. 


Figure  2.       Lattice  Prediction  Error  Filter. 


III.      THE  SPLIT-LEVINSON  ALGORITHMS 

The  split- Levinson  algorithms  are  based  on  the  theory  of  symmetric  and  antisym- 
metric polynomials.  We  know  that  for  an  k-th  order  real  autoregressive  process,  the 
normal  equations  are 


Rx     Rq 

R,       R^ 


R, 


•    Rk 

l 

a 

■   Rk-\ 

a\ 

0 

■   Rk-2 

a2 

= 

0 

■    R* 

<*k 

0 

(3.1) 


Rak  =  [a,0,0,.  ...0]7" 


(3.2) 


Using  the  Barlett  Bisection  Theorem  [Ref.  5:  pp.  1074-1076],  and  because  of  the  sym- 
metry of  the  autocorrelation  matrix,  we  can  say  that  the  predictor  coefficient  vector  is 
the  linear  combination  of  a  symmetric  and  antisymmetric  predictor  vector  given  by 


a    -a(5)  +  a(a) 
*k  -  zk    +  *k 


(3.3) 


The  symmetric  and  antisymmetric  vectors  are  defined  as 


a(5)- 
a  k    - 

'  B~ 

_#'_ 

a(a)- 

"   B 
-B' 

(3.4) 


where  B  represents  one-half  of  the  vector  components  of  a*,  and  B'  represents  the  re- 
versal of  the  vector  components  of  B.  Using  Eqs.  (3.3)  and  (3.4),  we  can  transform  the 
normal  equations  into 


Raf  =  CW2.0.0 W2]r 

*if>-  0/2,0,0 -<r/2]r 

Therefore,  we  can  see  some  favorable  consequences  of  these  revised  normal 
equations,  and  their  solutions.  First,  either  the  symmetric  or  antisymmetric  form  will 
give  the  same  solution,  and  second,  because  of  the  symmetry  of  the  predictor  vectors, 
we  need  only  solve  for  one-half  of  the  predictor  coefficients. 

Similar  to  the  Levinson  algorithm  we  now  proceed  to  develop  the  split-Levinson 
algorithms  from  the  input  sequence  autocorrelation  function  and  the  predictor  error 
vectors. 

A.     SPLIT-LEVINSON  ALGORITHM 

The  predictor  polynomial  ak(z)  is  defined  as 


«*(*)"/««*  '  (3-6) 


relative  to  the  given  Toeplitz  matrix  of  autocorrelation  lags.  Denote  the  reverse  of  our 
predictor  polynomial  as  a(z)  =  z-*^(z_1),  and  the  predictor  polynomial  has  been  shown 
to  obey  the  recursion  [Ref  3:  pp.  156-157] 

ak(z)  =  ak_{{z)  +  pkz~la(z)  (3.7) 

and  the  reverse  polynomial  of  Eq.  (3.7)  is 

ak(z)  =  z~\k_,{z)  +  pkak_x{z)  (3.8) 

We  now  want  to  form  a  new  polynomial  from  the  given  predictor  polynomial  that 
will  form  the  basis  of  the  split-Levinson  algorithm.  It  is  desired  to  show  that  the  deter- 
mination of  the  coefficients  of  this  polynomial  will  allow  us  to  recover  the  original  pre- 
dictor polynomial,  and  at  the  same  time  be  more  computationally  efficient.  We  define 
the  symmetric  polynomial  as  Pk{z),  and  the  antisymmetric  polynomial  as  P[c){z)  ,  and  we 
desire  them  to  be  of  the  form  [Ref.  1:  p.  472] 
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ft 
Pk(z)  =  £  PkiZ~' 

(3.9) 

;=0 

Recall  from  Eq.  (3.4)  and  (3.5)  that  the  symmetric  and  antisymmetric  predictor  coeffi- 
cients are  composed  of  two  vectors  that  are  reverses  of  each  other,  and  we  will  define 
these  vectors  so  that  they  obey  the  relationships 

Pk.i  =  -  Pkjc-l 

Consider  the  mathematical  interpretation  of  making  the  autocorrelation  matrix,  Rk 
a  singular  matrix.  If  the  reflection  coefficient  pk  is  made  ±  1,  then  this  corresponds  to 
an  element  of  Rk  making  the  matrix  singular.  For  this  reason  we  shall  designate  the 
symmetric  and  antisymmetric  predictor  polynomials  as  singular  predictor  polynomials 
[Ref.  2:  p.  472]  and  from  Eq.  (3.9)  they  are  defined  as 

Pk(z)  =  ak_,(z)  +  z-lak^(z) 

P{k%)  =  ak_l(z)-z-lak_](z) 

Also,  these  singular  predictor  polynomials  are  self-reciprocal  [Ref.  2:  p.  472]  because  of 
their  symmetry  and  may  be  expressed  in  the  following  forms 

Pt(,)=z-^(z-') 

From  Eq.  (3.11)  we  have 

z-'ak_x{z)  =  Pk{z)-ak_x(z) 
=  a,_,(z)-/f)(z) 

If  we  add  Eqs.  (3.7)  and  (3.8),  and  make  a  substitution  from  Eq.  (3.13)  we  have 

ak(z)  +  ak(z)  =  z~^ak_x{z)  +  pkak_{(z)  +  ak_x(z)  +  Pkak_x{z) 

=  (1  +  Pk)ak_,{z)  +  (1  +  Pk)2-'ak_x{z)  (3-14) 

=  >-kPfc) 
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where  we  have  defined  ).k  as 

;-*-1+Pa  (3-15) 

In  a  similar  fashion  we  can  solve  for  the  antisymmetric  normalized  singular  predictor 
polynomial  by  subtracting  Eqs.  (3.7)  and  (3.8),  and  substituting  from  Eq.  (3.13)  we  have 

ak(z)-ak(z)  =  Aka)P{ka)(z)  (3.16) 

where 

4a)=l-P/<  (3.17) 

Similar  to  the  predictor  polynomial  ak(z),  we  can  define  the  singular  predictor  coef- 
ficient vectors  for  Eq.  (3.9)  as  [Ref.   2  :  p.  472] 


Vk  =  [Pko>Pki>  •••  >Pkk] 


(3.18) 


Since  we  want  the  split-  Levinson  normal  equations  to  be  of  the  form 

Rk*k  =  tokQ,0,...,0f  (3.19) 

KAaA  =  [0 A-,  okf  (3.20) 

then  from  Eq.  (3.14)  and  (3.16)  the  singular  predictor  polynomials  are 


akU)       ak(z) 


«),.  _  «k(z)        ak(-) 

'  4 


(3-21) 


Pi%) 


Since  ak(z)  is  a  polynomial  formed  from  the  predictor  coefficient  vector  that  is  a  solution 
to  Eq.  (2.3),  it  follows  that  Pk(z)  and  P\a)(z)  are  solutions  to  the  Toeplitz  system  described 
by  Eqs.  (3.19)  and  (3.20).  [Ref.  1:  p.  472] 
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**P*  -  ^A 


Rk?f  =  Rk 


(3.22) 


Normalizing  Eqs.  (3.19)  and  (3.20)  by  ).k  and  /.k'\  the  split-Levinson  normal  equations  in 
matrix  form  are  expressed  as 


^P^UAO,...,!,]7" 

KkPk  -  Lr/<  -0'0'  •••  '-T*  J 

where  we  have  defined  the  modified  prediction  error  norm  [Ref.  2:  p.  651] 


(3.23) 


°k 

Ak 

ja)  _  _Zk_ 

Tk    -    ,(a) 


(3.24) 


If  we  expand  the  matrix  expressions  in  Eq.  (3.23),  the  modified  error  norms  may  be  ex- 
pressed as  finite  sums  of  the  predictor  coefficients 


I* 


-  /_^   LH  Pki 

i=Q 
k 


(3.25) 


M) 


i* 


Pki 


where  #,  is  the  i-th  autocorrelation  element  of  the  k-th  sub  matrix. 

Since  the  symmetric  and  antisymmetric  polynomials  are  closely  related,  we  shall 
derive  only  the  symmetric  polynomial  recursion  equations,  and  then  simply  present  the 
results  for  the  antisymmetric  case. 

The  final  step  in  the  derivation  is  to  derive  a  three  term  recursion  formula  for  the 
symmetric  polynomial.  From  Eq.  (3.11)  and  (3.14)  we  have  the  surprising  result  that  the 
predictor  polynomial  ak(z)  can  be  obtained  from  a  linear  combination  of  successive  sin- 
gular predictor  polynomials  [Ref.  2:  p.  472].  First,  form  Pk^{z)  from  Eq.  (3.11),  and  then 
eliminate  ak(z)  using  Eq.  (3.14) 
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=  ak(z)  +  z-][s.kPk(z)-ak(z)] 

-i  i  (3.26) 

=  (l-z   [)ak(z)  +  z   lAkPk(z) 

(\-z-x)ak(z)  =  Pk+l(z)-z-';.kPk(z) 

If  we  replace  k  by  k-1  in  Eq.  (3.26)  above  we  also  have 

(1  -  z-l)a*-iM  =  P&)  ~  z~\-iPk-i(*)  (3-27) 

We  now  form  our  recursion  formula  by  mutiplying  Eq.  (3.1 1)  by  (1  -  z~l),  and  use  Eqs. 
(3.15),  (3.26),  and  (3.27)  to  eliminate  pk  and  all  ak  predictor  polynomials. 

(1  -  z-l)ak(z)  =  (1  -  z-[)ak_,(z)  +  pk{\  -  z"V  Vito 

=  (1  -  z-l)ak_{(z)  +  p,(l  -  z-l)lPk(z)  -  ak_x{z)-] 

=  (1  -  z-l)ak_,(z)  +  pk[Pk(z)  -  ak_x{z)  -  z~lPk(z)  +  z_  Vito]  (3-28) 

=  «k-\W  -  z~l  ~Pk  +  ^Pk] 

=  a,_1(z)[(l-z-1)(l-p,)] 

If  we  now  substitute  for  (1  -V'K(^').  (1  -  z^K-ite-1)  from  Eqs.  (3.27)  and  (3.28),  and 
eliminate  pk  using  Eq.  (3.15),  we  can  complete  the  derivation 

Pk+l(z)  -  z~lAkPk(z)  =  lPk(z)  -  z-l;.,_,JP,_,(z)](2  -  Ak)  +  [AkPk(z)  -  Pk(zW  -  z"1) 
=  2Pk(z)  -  ).k{z)Pk{z)  -  2z"1;.,_lJP,_1(z)  +  z-,/,;.,_1P,_,(z) 
+  AkPk(z)  -  Pk(z)  -  z~]AkPk(z)  +  z~lPk(z)  (3.29) 

Pk+](z)  =  (1  -  z-y,(z)  4-  z"  VA-i(z)C4  "  2] 
=  (l+z-y,(z)-^z-1P,_1(z) 

Taking  the  inverse  Z  transform  of  Eq.  (3.29),  we  have  the  three  term  recursion  formula 
for  the  singular  predictor  coefficients 

Pki  =  Pki  +  Pkji-i  ~  *kPk-\ji-i  (3-30) 

where  the  recursion  parameter  <xk  is  defined  as 

Ofc-JjwB-iJ  (3.31) 

We  note  that  t*  is  determined  from  Pk(z)  from  Eq.  (3.23),  and  therefore  we  conclude  that 
the  singular  predictor  polynomials  can  be  recursively  computed  from  Eq.  (3.30).  How- 
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Poo  = 

0 

P\0  ~ 

1 

P\\     = 

-1 

*0 

ever,  the  recursion  parameter  ak  is  not  quite  in  the  correct  form.     From  Eqs.  (2.10), 
(3.15),  and  (3.31)  we  can  alternatively  compute  a*  from  [Ref.  2:  p.  473] 

*k—£?  (3.32) 

The  dual  relationships  for  the  antisymmetric  split- Levinson  formulas  can  be  derived 
by  following  a  procedure  similar  to  the  one  presented  above.  It  suffices  to  replace  the 
quantities  ).k%  xk<xk, pki,  by  their  antisymmetric  duals,  i.e.,  ptf,  and  use  the  following  anti- 
symmetric initial  conditions.  [Ref.  2:  p.  649] 


(3.33) 


Recursive  equations  for  the  symmetric  split-Levinson  -algorithm  are  summarized  in 
Table  4  of  Appendix  A.  Examining  the  entries  in  Table  4,  we  see  that  a  full  iteration 
loop  of  the  algorithm  requires  approximately  t  real  multiplications.  However,  because 
of  the  symmetry  of  the  singular  predictor  coefficients,  we  only  have  to  perform  one-half 
of  these  calculations.  Therefore,  for  a  k-th  order  filter  we  need  to  make  on  the  order  of 
k2jl  real  multiplications.  The  S  function  in  Table  4  is  used  to  distinguish  between  even 
and  odd  orders  of  the  indexing  variable. 

The  FORTRAN  program  SPLIT,  in  Appendix  B,  estimates  the  predictor  coefficients 
using  the  Levinson  and  split-Levinson  algorithms.  Figure  3  is  a  graphical  comparison 
between  the  known  test  filter  coefficients  of  SPLIT,  shown  by  the  solid  curve,  and  the 
filter  coefficients  computed  by  the  Levinson  and  split-Levinson  algorithms,  shown  by  the 
dashed  curve.  We  now  undertake  the  derivation  of  the  lattice  form  of  the  split-Levinson 
formulas  to  verify  the  numerical  complexity  of  that  method,  and  to  investigate  the 
symmetric  and  antisymmetric  lattice  structures  compared  to  the  Levinson  lattice  forms. 

B.      SPLIT  LATTICE  ALGORITHM 

We  begin  the  split  lattice  derivation  by  introducing  the  symmetric  and  antisymmetric 
error  vectors  xk,  \\c)  [Ref.  2:  p.  648].  If  we  use  previously  established  singularity  concepts, 
and  substitute  ±  I  for  pk  in  Eq.  (2.20)  for  the  symmetric  and  antisymmetric  error  vectors, 
x*  and  xS?\  respectively,  we  have 
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Figure  3.       Levinson  /  Split  Levinson  Coefficient  Comparison. 
**W-/*-i(')  +  *a-i('-1) 


(3.34) 


As  in  the  split- Levinson  case,  we  shall  proceed  with  the  derivation  of  the  symmetric 

split  lattice,  and  present  the  antisymmetric  lattice  results  at  the  end  of  the  derivation 
with  any  significant  changes  noted. 
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If  we  extend  the  singular  polynomial  concept  to  the  singular  predictor  coefficients, 
we  can  start  with  the  Levinson  coefficient  recursion  formula  and  substitute  ±  1  for  pk 

aki  =  ak-u  +  Pk.ak-\,k-i  (3-35) 

and  substituting  for  pk 

Pki  =  ak-l,i  +  ak-ljc-l  (3'36) 

If  we  write  Eq.  (3.34)  representing  the  time  index  (t)  with  the  subscript  (i)  we  have  an 
algorithm  that  is  more  easily  adapted  for  computers. 

xw/k-ut  +  bk-Ht-i  (3-37) 

Now,  comparing  Eqs.  (3.36)  and  (3.37)  we  have  a  direct  correlation  between  the  two, 
and  from  the  split-Levinson  equation  for  the  forward  error  vector,  we  can  write  the  dual 
split  lattice  equation  for  xK(t) 


(r)  =  £^s(/-0  (3.38) 


xk 

1=0 

Since  Eq.  (3.38)  is  in  the  form  of  a  convolution  sum  we  can  apply  Z  transform  theory 
to  see  if  any  inferences  can  be  made 

Xk(z)  =  Pk(z)S(z) 

S(z)     ~Pk{2) 

From  Eq.  (3.39)  we  can  conclude  that  the  symmetric  polynomial  constitutes  the  Z 
transform  of  the  transfer  function  formed  from  the  error  vector  and  input  sequence  z 
polynomials.  Now  we  can  use  the  previously  derived  split-Levinson  algorithm,  and  in- 
verse transform  it  to  obtain  the  lattice  error  vector  recursion  algorithm. 
Repeating  the  split-Levinson  recursion  we  have 

Pk+l(z)  =  (1  +  z-])Pk(z)  -  *kz-'Pk_x{z) 

=  Pk(z)  +  z-lPk(z)-cckz-]Pk_l(z) 

Multiplying  Eq.  (3.40)  by  S(z)  and  using  Eq.  (3.39)  for  Pk{z)S{z)  we  have 
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S(z)Pk+l(z)  =  S(z)lPk{z)  +  z~lPk(z)  -  xkPk_x(z)l 
Xk+l(z)  =  Xk(z)  +  z~]Xk(z)  -  zkz-xXk_x{z) 

Applying  the  inverse  Z  transform  to  each  side  of  Eq.  (3.41)  we  have  the  singular  pre- 
dictor error  vector  recursion  formula  [Ref.  2:    p.  650] 

W)  =  **(0  +  -"'^(0  -  <*fc*" Vi(')  (3-42) 

The  symmetric  lattice  structure  given  by  Eq.  (3.42)  is  shown  by  Figure  4.  [Ref.  2: 
p.  650] 

From  Eq.  (3.32)  we  know  that  the  recursion  parameter  ak  is  defined  as  rkjxk_u  and 
since  a.k  appears  in  the  recursion  formula  for  the  singular  error  vector,  we  need  to  solve 
for  it.  We  begin  with  the  Levinson  error  norm  equation 

*fc=(l  -Pk)<>k-\ 

2ak  =  2ak_x{\  +  pk)(\  -  pk) 
,      °k  ,        n         ,  (3.43) 

2ak_{{l  -  pk)  =  2rk 

where  we  have  substituted  rk  from  Eq.  (3.32),  and  2crJ(r_l(  1  -  pk)  is  defined  as  ||jcJ|2  [Ref. 
2:  p.  650]. 

The  initial  conditions  for  Eq.  (3.42)  must  be  examined  because  most  cases  are  trivial 
except  for  the  case  of  k.  =  0.  From  Eq.  (3.38)  we  have 

x0(t)=Poos(t),  (3.44) 

ind  from  [Ref.  2:  p.  648]  we  define  p00  =  2. 

The  antisymmetric  duals  are  very  similar  to  the  symmetric  case,  and  can  be  formed 
by  replacing  the  symmetric  variable  by  its  antisymmetric  counterpart,  i.e.,  x\a)(t)  for 
xk(i),  -  pk  for  pk,  etc.  The  initial  conditions  for  the  antisymmetric  case  are  given  below 

i»s  (3-45) 

xka\t)  =  s{t)-s{t-\) 

The  antisymmetric  lattice  structure  given  by  the  antisymmetric  dual  of  Eq.  (3.42)  is 
shown  in  Figure  5.  The  split-Levinson  lattice  formulas  given  by  Eqs.  (3.37),  (3.43),  and 
(3.44)  are  summarized  in  Table  5  of  Appendix  A.  If  we  examine  the  number  of  real 
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Figure  5.       Antisymmetric  Lattice  Structure. 
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multiplications  given  in  Table  4  and  Table  5,  then  we  can  deduce  the  following  conclu- 
sions. The  split  versions  of  the  Levinson  algorithm  do  produce  a  reduction  in  the  com- 
plexity of  the  calculations  when  compared  to  the  classical  versions.  The  split-Levinson 
produces  a  reduction  of  one-half,  and  the  Levinson  produces  a  reduction  of  3/2  [Ref.  2: 
p.  645].  The  FORTRAN  program  SLATIS,  Appendix  C,  implements  the  Levinson  and 
split-Levinson  lattice  algorithms,  and  a  graphical  comparison  between  the  known  test 
coefficients,  shown  by  the  solid  curve,  the  coefficients  estimated  by  the  Levinson  lattice 
algorithm,  shown  by  the  dashed  curve,  and  the  coefficients  estimated  by  the  symmetric 
split-Levinson  algorithm,  shown  by  the  dotted  curve,  is  presented  in  Figure  6. 

The  split  lattice  structures  shown  in  Figure  4  on  page  19,  and  in  Figure  5  on  page 
20  show  that  the  classical  lattice  structure  appears  to  be  lost  in  the  new  algorithm.  The 
distinct  advantage  of  the  original  lattice  structure  is  the  modularity  of  the  filter.  In  order 
to  retrieve  this  appealing  feature  of  the  lattice  filter  we  shall  now  proceed  to  derive  a 
revised  version  of  the  split  lattice  structure,  and  see  if  it  can  have  a  form  similar  to  the 
classical  structure. 

C.     SPLIT  LATTICE  REVISED  STRUCTURE 

To  begin,  consider  the  second  order  classical  lattice  structure  derived  from 
Figure  2  on  page  8.  We  can  write  the  following  equations  for  the  first  and  second  stage 
forward  and  backward  prediction  errors, 

/1(n)  =  s(«)  +  p,j(«-l) 
S\{n)  =  pAn)  +  s(n  -  1) 

Mn)=Un)  +  p2g[(n-\)  (J>Ab) 

gi{n)  =  P/i(«)  +  g\{n  -  1) 

Now  substituting  the  equations  for/(n)  and  forg,(n)  into  the  equations  for/2(n)  and 
g2(n).  solving  for  the  second  stage  forward  and  backward  prediction  filter  errors  and 
taking  the  Z  transforms  yields  the  transfer  functions 

-^f-  =  1  +  z_1(p,  +  P]p2)  +  p2z~2  (3.47) 


and. 


-|jg^  =  p2  +  z~\px  +  plPl)  +  z~2  (3.48) 
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Figure  6.       Levinson  vs.  split-Levinson  Coefficient  Comparison. 

which  obviously  yield  the  second  order  predictor  transfer  function  A2(z)  and  its  reverse 
version  A2(z),  where  a20  =  1,  an  =  p2.  a^  =  p,  -  pxp2  =  p,(l  -  p2).  Now  forming  the  third 
order  symmetric  polynomial  P%{z)  from  our  second  order  example,  we  have  [Ref.  1:  p. 
472] 


P3(:)  =  A2(z)  +  z-2A2(z~') 


ay  +  (a,  +  a2)z   '  +  (a,  +  a2)z  2  +  %z 


(3.49) 


If  we  define  a0  =    1.  and  («,  4-  a2)  =  /?,  .  then 


P3(r)-l+/>1r-,+/.1z-2  +  2-J 


(3.50) 


Define  the  followinc  terms: 


F1(z)=l+Pl2-' 
G1U)  =  r-,(l  4-/7,2) 


(3.51) 


therefore, 


^(r)  =  JF,(r)  +  z-2G,(r) 


(3.52) 


Equation  (3.52)  defines  the  revised  symmetric  split-lattice  structure,  and  Figure  7  gives 
a  graphical  representation  of  that  structure. 


Figure  7.        Proposed  Split  Lattice  Configuration 


In  order  to  show  that  the  preceding  classical  structure  can  be  equated  with  the 
symmetric  polynomial  derived  earlier,  it  is  now  necessary  to  form  the  forward  and 
backward  prediction  error  transfer  function  of  the  new  split  lattice  structure  and  com- 
pare them  to  the  Z  transform  of  the  symmetric  polynomial.  Therefore,  from  Figure  7 
and  Eq.  (3.52),  we  can  write  P3(z)  as 


P3{z)  =  1  +  A',2"'  +  z"2(A',  +  z"1) 


(3.53) 
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Now,  if  we  compare  Eq.  (3.50)  to  Eq.(3.53),  we  see  that  K}  -  pv  But,  from  Eq.  (3.49) 
we  know  that  px  =  (a,  +  a2),  and  substituting  for  a{  and  a2  from  Eq.  (3.48)  and  (3.49) 
yields 

Kl=px=px+p1(l-pl)  (3.54) 

Let  us  now  rederive  the  symmetric  polynomial  from  the  revised  split  lattice  structure. 
From  Eq.  (3.50)  and  (3.54)  we  have 

P3{z)  =  1  +  (p,  -  p2  -  p,p2)z~X  +  z~2[(p,  -  p2  -  pxPl)  +  z_1]  (3.55) 

Since  we  have  found  that  the  symmetric  lattice  can  be  restructured  to  a  form  similar 
to  the  classical  lattice,  the  next  logical  step  is  to  find  a  recurrence  relation  for  the  new 
lattice.  Let  us  consider  a  5;/'  order  symmetric  polynomial  to  determine  the  step-down 
recursion  procedure,  given  by 

P5(z)  =  l+p5,z-X  +  p52z~2  +P5iz~3  +Ps^~4  +  z~5  (3-56) 


From  Eq.  (3.51)  we  have 


F2(z)=\+P5,z   l+p52z  2 
G2(z)  =  z"2+/>51z"'  +p52 


(3.57) 


As  per  the  observations  made  in  Figure  7  on  page  23  and  Eq.  (3.51),  we  have  the  lattice 
reflection  coefficient  for  the  second  stage,  K2  =  pS2.  Now  reduce  the  order  of  F(z)  to  find 
the  First  stage  reflection  coefficient  using  the  standard  inverse  Levinson  recursion  [Ref. 
4:  pp.  156-157]. 

F2(z)-K2G2(z) 

i(z)= — r~j? — 

1  "  Kl  (3.58) 

=  1+     P5]      -- 


therefore 


1  +  A', 


k,=ttV  (3-59) 


Now  rewriting  the  equations  for  Fx{z)  and  F2{z)  using  the  derived  reflection  coefficients, 
we  have 
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and  we  know  that 


f,U)=l  +  A'1z-1 

Gl{:)  =  z-1  +  A', 

F2(z)  =  1  +  A*,(  1  +  K2)z~]  +  K2z'2 

G2(z)  =  z~2  +  A,(l  +  AOr"1  +  K7 


P5{z)  =  F2(z)  +  z-3G2(:) 


(3.60) 


=  1  +  A,  ( 1  +  K2)z~]  +  A2z~2  +  A2z"3  (3.61) 

+  A',(l  +  A2)z~4  +  z~5 

Equating  the  terms  in  Eq.  (3.56)  and  (3.61),  we  have  p5l  =  AT,  (I  +  A,)  and  p<2=  K2. 
Knowing  the  values  of  A;  and  K2,  we  now  can  form  a  two  stage  symmetric  lattice  similar 
to  that  in  Figure  7  on  page  23  to  implement  Ps{z)  .  However,  we  are  interested  to  find 
if  we  can  recursively  obtain  the  lower  orders  P4(z),  Pj(z),  etc.  or  the  higher  orders 
P6(z),  P-(z),  etc.  from  Ps  (z).  In  an  attempt  to  form  P6(z),  we  use  the  standard  forward 
recursion  [Ref.  3:  pp.  156-157]  to  obtain 

F3(r)  =  f2(;)  +  A3z-1G2(z) 

=  1  +  A',(l  +  A"2)z_1  +  K2z~2  +  K2K3z~l  (3.62) 

+  AiA'3(l  +  K2)z~2  +  K3z~3 


and 


g3(z)  =  z-3  +  ^(i  +  a:2)2-2 


+  AT2z~'  +  A2A3z"2  (3.63) 

+  K{K3{1  +  K2)z~l  +  K3 

Forming  the  symmetric  polynomial  P6(z)  from  Eqs.  (3.62)  and  (3.63)  we  have 

P6(z)  =  F3(z)  +  z-3G3(z) 

=  1+(A,  +  A,A2  + A2A3)z_I 

+  (K2  +  A,A"3  +  KxK2K3)z-2  +  K3z~3  (3.64) 

+  A3z"3  +  (K2  +  A,  A3  +  A,A2A3)z"4 

+  (K]  +  K:K2  +  K2K3)z~5  +  z~6 

Comparing  Eq.  (3.64)  to  the  symmetric  form  of  P6(z),  we  have 


p6l  =(A',  +  KXK2  +  K2K2) 

^51  /,5l/,52      ,     1  (3.65) 

The  one-half  enters  into  Eq.  (3.65)  because  we  know  that  for  even  orders  the  polynomial 
coefficients  are  symmetric  about  the  center  element,  and  they  must  be  shared  in  the 
matrix  equations.  We  shall  now  expand  Eq.  (3.65)  to  attempt  to  develop  a  recursive 
algorithm  for  the  sixth  order  coefficients  from  the  fifth  order  coefficients.  Expanding  Eq. 
(3.65)  we  have 

/  1  +  p52  \ 

0  +Psi)P6\  =  Ps\  +Ps\P52  +  [  2 )&n  +  Psi  ~  "sPul  (3.66) 

where  the  term  in  brackets  is  an  expansion  of  the  coefficient  recursion  formula,  Eq. 
(3.30),  for  p6i.  Collecting  terms  we  have 

(1  +P52)P6\=Ps\^+P52)  +  ^+P52)T^lP52+P52-^5P^'] 

oc5  (3-67) 

Pb\  =  P5\  +  P52\j>52  ~  ~YP42~1 

Substituting  for/?J2  from  Eq.  (3.30) 

PS2=P42+P*\  ~a<lP3\  (3-68) 

From  Eq.  (3.68)  we  can  observe  that  the  ak  recursion  parameter  and  the  number 
of  previous  coefficients  required  to  be  known  are  increasing  in  order,  and  it  appears  that 
a  simple  recursive  algorithm  based  on  the  above  approach  is  not  possible.  Note  that 
although  the  new  lattice  structure  does  not  appear  to  be  order  recursive,  we  can  express 
a  given  order  symmetric  or  antisymmetric  lattice  structure  in  a  more  conventional  form. 
In  summary,  we  know  that  the  split-Levinson  algorithm  is  a  viable  replacement  for 
the  classical  algorithm  because  of  its  computational  efficiency.  We  have  studied  both 
autocorrelation  and  data  (or  lattice)  realizations  of  the  split-Levinson  algorithm.  An  at- 
tempt to  derive  a  recursive  split  lattice  algorithm  yielded  a  classical-like  lattice  structure, 
but  it  is  not  recursive  in  order.  Further  investigation  is  necessary  in  this  direction.  We 
now  need  to  test  the  split-Levinson  algorithms  on  some  signal  processing  applications. 
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IV.      APPLICATIONS  OF  THE  SPLIT-LEVINSON  ALGORITHM 

In  this  chapter,  we  apply  the  split- Le Vinson  algorithm  in  (1)  the  MA  parameter  es- 
timation problem,  and  (2)  the  extended  Prony  method  of  spectral  line  estimation.  Before 
we  take  up  these  two  applications  we  examine  the  algorithm's  usefulness  if  the  original 
filter  has  coefficient  symmetry,  i.e.,  the  impulse  response  of  a  linear  phase  FIR  filter. 

A.      HANKEL  AND  TOEPLITZ  MATRICES 

In  previous  derivations  we  have  assumed  the  FIR  filter  equation  to  be  non- 
symmetric.  Let  us  now  investigate  the  problem  where  the  Filter  equation  is  symmetric, 
i.e..  of  the  form 

y{n)  =  s(n)  +  a{s(n  -  1)  +  a2s(n  -  2)  +  ••• 

(4.1) 
+  ak^s(l)  +  aks{n  -  k) 

By  definition,  a  symmetric  polynomial  is  self-reciprocal,  that  is 

ak(z)  =  ak(z)  =  z-kak(z-[)  (4.2) 

Therefore,  from  the  Levinson  algorithm,  predictor  polynomials  are  known  to  obey  the 
recurrence  relation 

ak(z)  =  ak_x{z)  +  pkz~]ak_}(z)  (4.3) 

and  in  our  special  case  we  have 

ak(z)  =  ak_{{z)  +  pkz^ak_{{z)  =  a(z) 

_i  (4-4) 

=  d  +  p,z  ViU) 

In  order  to  formulate  a  set  of  equations  similar  to  the  split-Levinson,  it  is  necessary 
to  derive  the  normal  equations  for  our  special  case,  and  compare  them  to  the  standard 
equations,  in  order  to  develop  the  recursive  algorithms.  Since  the  predictor  coefficients 
in  a  recursive  algorithm  produce  estimates  of  s(n)  as  the  algorithm  steps  through  its  re- 
cursive steps,  denote  this  estimate  as  s{n).  In  vector  form,  we  then  have 
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s{n\n  -  1)  =    -  [s(/7  -  1)  s{n  -  2)  ...  5(0)] 


«2 


(4.5) 


To  derive  the  normal  equations  we  must  find  the  minimum  mean  squared  error  (MSE) 
from  the  equation  for  the  error, 


e(n)  =  s(n)  —  s{n\n  —  1) 


(4.6) 


The  minimum  mean  squared  error  is  found  by  squaring  the  error  term,  and  then  differ- 
entiating the  squared  term  with  respect  to  the  given  ak  vector.  Combining  these  two  ev- 
olutions we  have  the  following  equations, 


MSE  =  J 

=  E\.e\n)-] 

=  £[(*(*) -(5>)|*-1))2] 

To  obtain  the  normal  equations,  we  carry  out  the  following  steps: 

jL  =  o  =  -2£[5(«)s[_,]  -2£Is[_,si_1] 


(4.7) 


(4.8) 


From  the  split- Levinson  recursion  formulas  we  know  that  the  singular  symmetric 
polynomial,  Pk{z),  is  defined  as  the  following, 


Pk{z)  =  ak_x(z)  +  z   V,(z) 


(4.9) 


Since  our  predictor  polynomial  is  symmetric,  it  is  a  reasonable  question  to  ask  if  sym- 
metric polynomials  also  obey  this  recursive  relationship.  It  is  noted  as  an  immediate 
consequence  of  the  recursive  problem,  all  of  the  preceding  polynomials  will  also  be 
symmetric.  Therefore,  we  check  the  singular  predictor  polynomial  recursion  to  see  if  it 
holds  when  the  original  polynomial  is  itself  symmetric 
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ak(z)  =  aA_j(z)  +  z  Xak_{{z) 
=  ak_x(z)  +  z-'z-\_x{z) 

=  1  4-  a,;'1  +  a2z~2  +  ...  +  z~{k~])  +  z~k[\  +  a,z~x  +  a2z~2  +  ...  +  z~k]    l        J 
=  1  +  (1  +  a,);"1  +  (a,  +  a,)z~2  +  ...  +  z~* 

Now  that  we  see  that  the  Levinson  recursive  equation  holds  for  a  symmetric 
polynomial,  we  derive  the  recursive  relationships  for  our  polynomial  using  what  is 
known  from  the  split- Levinson  equations.  We  have  defined  the  symmetric  polynomial 
Pk(z)  to  be  a  normalized  combination  of  a  non  symmetric  polynomial,  ak(z),  and  its  re- 
ciprocal image,  ak(z)  in  the  form  of, 

AkPk(z)  =  ak(z)  +  ak(z)  (4.11) 

By  direct  substitution  it  is  a  trivial  matter  to  show  that  this  relationship  also  holds  for 
a  symmetric  polynomial,  ak{z). 

In  order  to  develop  the  recursion  for  the  symmetric  polynomial,  it  is  necessary  to 
express  the  desired  linear  predictor  in  terms  of  the  previous  two  predictors.  To  this  end 
use  Eq.  (4.10),  to  form  ak^x{z)%  and  substitute  from  Eq.  (4.11)  to  perform  this  task. 

ak(z)  =  ak_x(z)  +  z~{  a^iz) 
ak+{(z)  =  ak(z)  +  z-lak(z) 

=  ak(z)  +  z~\/.kak(z)  -  ak{z)~] 
=  ak(z)(\-z-])  +  z-').kak(z) 

Solving  for  (1  —  z~l)ak(z),  and  forming  the  quantity  (1  —  z~l)ak_x(z)  we  have 

(1  -  z~X)ak{z)  =  ak+x(z)  -  z~X).kak{z) 
(1-z     K_,(z)  =  ak(z)  -  z     ).k_xak_x(z) 

These  relationships  will  now  allow  us  to  form  the  three  term  recursion  for  the  given 
symmetric  polynomial  from  Eqs.  (4.3),  (4.11),  (4.12),  and  (4.13) 

ak{z)  =  ak_{(z)  +  pkak_{(z) 

-i  _i  (4.14) 

ak+\(z)  =  0  +z     )ak(z)  ~  akz     ak-\(z) 

where  we  have  defined  ak 

«k  =  >.k_x{2-Xk)  (4.15) 
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From  Eq.  (4.14)  we  can  see  that  the  coefficient  recursion  formula  is  the  same  form 
as  Eq.  (3.29),  and  we  can  deduce  that  the  split- Levinson  algorithms  will  work  equally 
well  for  symmetric  polynomials  that  describe  unknown  filters  as  it  does  for  polynomials 
that  are  not  symmetric. 

B.     FIR  MOVING  AVERAGE  PARAMETER  ESTIMATION 

If  we  consider  an  FIR  filter  with  an  input  sequence  given  by 
sr=  [s(n)s{n  —  l)...s(n  —  m)],  and  an  output  y(n)  given  by 

M 

y(n)  =  Yjais(n-i)  (4.16) 

(=0 

then  we  can  develop  the  necessary  equations  to  estimate  the  moving  average  parameters, 
and  solve  for  the  FIR  filter  coefficients.  The  algorithm  to  estimate  the  predictor  coeffi- 
cients can  be  defined  as  follows: 

Let  the  three  predictions,  yf{n),  sf{n),  and  s^(n),  represent  the  m-th  order  predictions 
of  the  forward  estimate  of  y,  the  forward  estimate  of  s,  and  the  backward  estimate  of  s 
respectively.  [Ref.  6:  p.  1  ] 


#"(«)  =  Xm«-0  (4-17) 

i=0 


(«)  =  £C/5(«-0  (4.18) 


■y 

/=0 


$"(„  -m)  =  J^di  s(n-m  +  /),  (4. 19) 

where  the  coefficient  vectors  are  defined  as, 

br=[l-*o  -b,...-bm] 

cr=[0  1  -c,...-cj  (4.20) 

dr=[0-cm  -cm_,...-cx^ 

Without  going  through  all  the  details  for  the  MA  parmeter  recursions,  we  can 
summarize  the  recursion  formulas  for  the  predictor  coefficients  as  [Ref.  6:  p.  2] 
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h        \ 

0  J  +  W 


Notice  that  the  predictor  vector  dm  is  not  included  in  the  preceding  definitions  since  it  is 
the  reverse  c77 

If  we  examine  Eq.  (4.21)  we  see  that  the  recursive  relationship  for  cm  is  a  statement 
of  the  Levinson  recursion,  since  K"t  and  K™  are  the  m-th  order  reflection  coefficients. 
Therefore,  we  can  apply  the  split-Levinson  algorithm  to  solve  for  c"\  form  &"  ,  and 
recursively  determine  b"\  Finally,  from  the  theory  of  Moving  Average  processes,  b"7  =  zr 

The  FORTRAN  program  MAV1,  in  Appendix  D,  uses  a  25-th  order  FIR  test  filter 
to  to  generate  a  test  sequence,  and  the  results  are  given  in  Table  6  of  Appendix  A. 

Figure  S  is  a  graphical  comparison  between  the  known  test  coefficients,  shown  by 
the  solid  curve,  and  the  computed  filter  coefficients,  shown  by  the  dashed  curve.  The 
curves  are  fitted  to  the  linear  magnitudes  of  the  coefficients  by  interpolating  spline 
techniques. 

C.      EXTENDED  PRONY  METHOD 

The  estimation  of  the  existence  of  sinusoids  in  the  presence  of  noise  is  a  common 
occurrence  in  signal  processing  applications.  In  this  simulation,  we  will  show  that  the 
split-Levinson  algorithm  can  be  directly  implemented  in  the  process  to  determine  the 
approximate  frequencies.  The  noise  is  zero  mean,  unit  variance,  and  white  in  nature. 

In  this  application  of  the  split-Levinson  algorithm,  we  attempt  to  approximately  fit 
p  exponentials  to  a  data  set  of  N  samples.  We  have  the  constraint  that  N  >  2p,  and  a 
least  squares  estimation  procedure  is  used.  We  begin  by  defining  the  estimate  of  our  set 
of  data  samples.  [Ref.  7:  p.  1404] 


-s 


bmznm     n  =  Q,..,N-\  (4.22) 


where  b„  =  Am  exp(/0J/2,  and  zm  =  exp(J2nfmAt).  The  zm  's  are  roots  of  unit  modulus  with 
arbitrary  frequency  and  occur  in  complex  conjugate  pairs  as  long  as  fm  #  0  or  l/2Af 
Therefore,  in  order  to  solve  for  the  p  sinusoids,  we  must  solve  for  the  roots  of  the  fol- 
lowing polynomial.  [Ref.  7:  p.  1406] 
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Figure  8.       MA  Coefficient  Comparison. 


¥(z)  =  Y    akz2p-k  =  0 


(4.23) 


The  roots  can  be  of  unit  modulus,  and  occur  in  complex  conjugate  conjugate  pairs  if 
we  constrain  the  polynomial  coefficients  to  be  symmetric  about  the  center  element  of  the 
polynomial  [Ref.  7:  p.  1407].  This  is  an  exact  ocurrence  in  the  symmetric  and  anti  sym- 
metric polynomials  of  the  split-Levinson  algorithm,  as  long  as  the  order  of  the 
polynomial  is  even,  that  is 


ir=|^la1  ...  Vi  ~j~  \ 


(4.24) 


Note  that  the  last  element  of  the  coefficient  vector  is  ap,2  rather  than  ap  because  of  the 
symmetry  of  the  polynomial,  and  that  symmetric  polynomials  only  guarantee  that  if  a 
root  z,  occurs,  then  so  does  its  reciprocal  z;1  [Ref.  7:  p.  1407]. 

The  program  EPRONY1,  Appendix  E,  utilizes  the  split-Levinson  algorithm  to  ap- 
proximate the  p  order  sinusoids  to  the  given  data  set.  A  summary  of  the  simulation  cases 
studied  are  presented  in  Table  1,  and  the  graphical  results  follow. 

Table   1.      SUMMARY  OF  TEST  CASES 


Case 

Number 

Constants 

Variables 

1 

Npts 
Test  frequencies 

fs 

SNR(-10,  0,  10  dB) 

2 

Test  frequencies 

/, 

SNR 

NPTS 

3 

Test  frequencies 
fs 

SNR 

Npts 

Filter  Order 

4 

Test  frequencies 

/, 

Npts 

SNR 

Filter  Order 

SNR  (-10,0,10  dB) 

1.  Simulation  Parameter  Definitions 

All  simulation  cases  are  done  in  the  presence  of  white  noise,  and  a  minimum 
possible  separation  frequency  for  the  input  sinusoids  was  determined.  Ail  plots  are 
sinusoid  magnitude,  in  a  linear  scale,  versus  digital  radian  frequency,  6  for  0  <  6  <  n. 

2.  Simulation  Results 

We  begin  the  spectral  line  estimation  simulations  with  all  parameters  fixed,  and 
then  selectively  choose  a  parameter  to  vary  and  observe  the  effects  of  these  changing 
parameters.  We  begin  with  two  sinusoids  of  fx  =  50  Hz,  f2  =  75  Hz,  which  yield 
0j  =  1.396  radians,  and  d2  =  2.0944  radians,  respectively.  The  number  of  data  points 
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(NPTS)  is  set  at  1500,  and  the  filter  order  (M)  is  chosen  to  be  4  indicating  the  presence 
of  two  sinusoids.  Now,  we  chose  to  vary  the  signal-to-noise  ratio  (SNR)  for  10  dB.  0  dB, 
and  -10  dB,  shown  in  Figure  9. 

Figures  9  (a)  and  9  (b)  show  that  lowering  the  SNR  from  10  dB,  Figure  9  (a), 
to  0  dB.  Figure  9  (b),  causes  the  estimation  to  worsen,  and  both  indicate  low  frequency 
estimation  error.  From  both  of  these  cases  we  can  deduce  the  presence  of  two  sinusoids, 
but  in  Figure  9  (c)  it  appears  that  only  one  sinusoid  is  present,  and  the  low  SNR  has 
caused  spectral  estimation  to  fail.  The  conclusion  for  the  first  case  is  that,  the  better  the 
SNR,  the  better  the  spectral  estimate. 

For  the  second  case  we  selected  NPTS  as  the  variable  parameter,  held  the  filter 
order  and  sinusoid  frequencies  constant  as  before,  and  set  the  SNR  at  0  dB.  From  Fig- 
ures 10(a),  (b),  and  (c)  we  clearly  see  that  the  better  estimation  occurs  with  NPTS  = 
1000,  Figure  10(b),  because  of  the  equal  amplitude  and  accurate  frequency  estimation 
as  compared  to  Figures  10  (a)  and  10  (c).  All  plots  show  low  frequency  error,  but  also 
suggests  that  simulations  should  be  conducted  with  NPTS  set  to  1000-1500  points  for 
the  best  results.  This  case  provides  the  rationale  for  the  value  of  NPTS  for  all  other 
simulations. 

In  the  third  simulation  all  parameters  are  fixed  as  in  the  second  case,  and  M  is 
varied  for  4,8,  and  10.  From  Figures  11(a),  (b),  and  (c)  we  can  see  that  although  there 
are  only  two  sinusoids  present,  the  estimation  plots  show  M/2  sinusoids  for  all  values 
of  filter  order.  Each  plot  has  frequency  components  in  the  vicinity  of  the  actual  fre- 
quencies, but  they  also  give  false  indications  of  spectral  lines.  If  we  were  making  a  deci- 
sion on  the  number  of  frequencies  based  on  the  magnitude  plots  of  Figures  11(b)  and 
(c),  then  signifigant  errors  would  be  introduced.  In  this  case  it  is  obvious  that  unless  the 
exact  number  of  sinusoids  present  is  known,  then  the  estimation  technique  will  fail. 

In  the  final  simulation  we  introduce  two  additional  sinusoids,  and  examine  the 
effects  of  SNR  on  spectral  line  estimation.  The  constants  for  this  case  are/  =  35  Hz, 
f2  =  S5  Hz,/,  =  105  Hz,/4=  175  Hz,  NPTS  =  1500,  M  =  4,  and/ =  525  Hz.  Digital 
frequencies  are  0.419,  1.010,  1.496,  and  2.094  radians  per  second  respectively.  As  in  case 
1,  SNR  takes  on  the  values  of  10  dB,  0  dB,  and  -10  dB.  In  Figure  12(a),  where  the  SNR 
is  10  dB,  we  get  good  results  with  near  uniform  amplitude  estimation,  and  digital  fre- 
quencies that  are  close  for  all  frequencies.  However,  from  Figures  12  (b)  and  12  (c),  we 
can  see  that  the  spectral  line  for/j  is  missing,  and  an  errant  line  appears  at  approximately 
2.8  radians.  Both  of  the  figures  show  unequal  amplitude  indications,  but  3  of  the  4 
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spectral  lines  are  in  close  proximity  to  their  actual  values.  From  this  case  we  draw  the 
conclusion  that  the  spectral  line  estimation  performance  deteriorates  at  low  SNRs. 
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Since  the  overall  purpose  of  the  simulation  was  to  test  the  applicability  of  the 
split-Levinson  algorithm  to  the  test  cases,  then  it  has  been  shown  that  the  split- Levins  on 
will  produce  estimates  for  the  respective  cases.  However,  if  we  examine  the  accuracy  of 
the  low  signal-to-noise  ratio  cases,  we  see  that  the  new  algorithm  suffers  a  similar  fate 
as  the  classical  case  in  that  it  is  difficult  to  accurately  estimate  the  correct  sinusoidal 
frequencies  in  the  presence  of  noise. 

D.     CONCLUSIONS 

The  split-Levinson  algorithm  has  been  shown  to  be  computationally  more  efficient 
than  its  classical  counterpart.  We  can  say  that  the  application  of  the  split-Levinson  al- 
gorithms to  practical  applications  in  lieu  of  the  Levinson  algorithm  can  be  advanta- 
geous, and  the  computational  cost  can  be  reduced  significantly  for  large  order  systems. 
Additionally,  the  split-Levinson  algorithms  are  applicable  to  problems  where  we  want 
to  model  MA  parameters,  and  perform  spectral  estimation  using  the  Prony  method. 

We  note  the  following  restrictive  areas  for  the  new  algorithm  that  could  make  it 
unsuitable  for  certain  signal  processing  applications. 

1.  Non-recursive  split-lattice  algorithm. 

2.  Computational  accuracy  degradation  when  performing  spectral  estimation  in  low 
signal-to-noise  cases. 

3.  Complexity  of  symmetric  and  antisymmetric  lattice  structures. 

Since  we  have  shown  that  the  split-Levinson  algorithm  is  a  viable  substitute  for  the 
Levinson  recursion,  it  is  reasonable  to  consider  areas  of  this  topic  for  further  research. 
We  know  that  the  symmetric  lattice  structure  can  be  expressed  in  a  classical  form  for 
given  filter  orders,  however,  a  recursive  algorithm  for  this  new  structure  has  been  elusive. 
We  propose  that  the  existing  recursive  algorithm  for  the  singular  predictor  coefficients 
should  be  studied  to  see  if  a  redefinition  of  equation  parameters  can  extract  a  new  re- 
cursive algorithm  for  the  new  symmetric  lattice.  Additionally,  the  algorithm's  poor  per- 
formance in  low  signal-to-noise  ratio  test  cases  of  the  extended  Prony  method  is  similar 
to  the  performance  of  the  classical  algorithm.  Therefore,  techniques  to  improve  the 
classical  algorithm's  performance,  as  detailed  in  [Refs.  8,9],  may  be  investigated  for  ad- 
aptation to  the  split-Levinson  algorithm. 
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APPENDIX  A.     TABULAR  SUMMARY  OF  ALGORITHMS 

The  tables  given  in  this  Appendix  were  taken  from  [Ref.  2:  pp. 648-674],  and  are 
presented  for  the  convenience  of  the  reader. 


Table  2.      THE  LEVINSON  ALGORITHM 


Computation 

Add 

Mult 

an,0  =  1 .         a0  =  c0 

For  k  =  l,2,...,n 

k-\ 

'■k  =  ~2-j    Ck-iak-\,i 
i=0 

k-1 

k-1 

ak,0  -  1.        Pk  =  *-kl°k-\ 

1 

ok  =  ak_x  -  Akpk 

1 

1 

akJ  =  ^k-\,i  +  Pkak-\Jc-i 

k-1 

k-1 

0=1,2 ,k-  1) 

41) 


Table  3.      THE  LEVINSON  LATTICE  ALGORITHM 

Computation 

Add 

Mult 

/0(r)  =  b0(t)  =  s(t)    [0<(<  X-  1) 

V-l 

a0  =  Y  s(t)2 

N-l 

N 

For  k=  1.2 n, 

X+k-2 

;=1 

N  +  k-3 

N  +  k-2 

PA  =  hl°k-\ 

1 

°k  =  °k-\  -  J-kPk 

1 

1 

A(0=/*-i(0  +  pA-i('-1) 

N+k-2 

N  +  k-l 

**(')  =  p&4-i(')  +  W) 

N  +  k-2 

N  +  k-l 

(t  =  0.1 N  +  k-1) 

Table  4.      THE  SPLIT-LEVINSON  ALGORITHM 


Computation 

Add 

Mult 

Po.o  =  2,  p,,0  =  pu  =  1,  t0  =  \  p0  =  0 

Fork=l,2,...,n 

fc+l  =  2r-<5, 

with  (5  =  0  or  1 

r-l 

(=0 

t 

/-  1  +3 

P*.o  =  !»   *k  =  TklTk-\ 

1 

Pk=  1  -  <**/(!+  Pi-i) 

2 

1 

Pft+U  =  fi«  +  fttf-l  ~  <*kPk-\J-\ 

2t 

t-1 

(/=  1,2,....,/) 
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Table  5.      THE  SPLIT  LATTICE  ALGORITHM 


Computation 

Add 

Mult 

x0{t)  =  2s(t)    (0<t<X-  I) 

x^t)  =  s(t)  +  s(t -  1)    (0<t<.\-) 

N-l 

,V-1 
P0  =  0  ,      T0  =  2]  5(/)2 
;=0 

N-l 

N 

For  k=  l,2,...,n 

,V+fc-l 

r=0 

N  +  k-1 

N  +  k 

a£  =  Tklrk-\ 

1 

Pi<=  1  -«*/(!+ P*-i) 

2 

1 

**+i(f)  =  **(')  +  xk(t  -  1)  -  a^_,(/) 

2(N+k-l) 

N  +  k-1 

(;  =  0, ...  ,iV+A) 
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Table  6.       MOVING  AVERAGE  TEST  RESULTS 

Coefficient 

Test 

Symmetric 

Index 

CoefYicients 

split-  Levinson 
Coefficients 

0 

0.1 

0.21 

1 

0.2 

0.12 

2 

0.3 

0.22 

3 

0.4 

0.32 

4 

0.5 

0.42 

5 

0.6 

0.52 

6 

0.7 

0.61 

7 

0.8 

0.71 

8 

0.9 

0.S1 

9 

1.0 

0.91 

10 

1.1 

1.02 

11 

1.2 

1.13 

12 

1.3 

1.23 

13 

1.4 

1.33 

14 

1.3 

1.24 

15 

1.2 

1.15 

16 

1.1 

1.06 

17 

1.0 

0.97 

18 

0.9 

0.88 

19 

0.8 

0.78 

20 

0.7 

0.69 

21 

0.6 

0.59 

22 

0.5 

0.50 

23 

0.4 

0.40 

24 

0.3 

0.31 

25 

0.2 

0.20 

26 

0.1 

0.10 
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APPENDIX  B.      SPLIT-LEV  INSON  PROGRAMS 


PROGRAM  TO  CALCULATE  THE  NTH  ORDER  FIR  PREDICTOR  FILTER  USING 
THE  SYMMETRIC  AND  ASYMMETRIC  SINGULAR  PREDICTOR  POLYNOMIALS, 
THE  SPLIT- LEVINSON  RECURSION  FORMULA,  AND  THE  AUTOCORRELATION 
FUNCTION  OF  THE  INPUT  SEQUENCE. 

VARIABLE  DEFINITIONS 

-  N-TH  DEGREE  NORM  OF  THE  FILTER. 

-  INTEGER  VARIABLE  USED  TO  CONTROL  ACCESS  TO 
SUBROUTINE  EVEN  WHEN  THE  INDEX  VARIABLE  K  IS  AN 
EVEN  INTEGER. 

-  REAL  VECTOR  USED  WHEN  DEFINING  THE  SINGULAR 
PREDICTOR  POLYNOMIAL  (PK(Z))  IN  TERMS  OF  THE 
NORMALIZED  SYMMETRIC  AND  ANTISYMMETRIC  PARTS  OF 
THE  PREDICTOR  POLYNOMIAL  AK(Z). 

LAMDA(K)  =  1  +  RHO(K) 
LAMDAS(K)  =  1.  -  RHOS(K) 

-  REAL  VECTOR  USED  TO  SIMPLIFY  THE  THREE -TERM 
RECURRENCE  RELATION  FOR  THE  SINGULAR  PREDICTOR 
POLYNOMIALS  OF  THE  SAME  TYPE. 

ALPHA(K)  =  LAMDA(K- 1)*(2.  -  LAMDA(K)),  OR 
ALPHA(K)  =  TAU(K)/TAU(K-1) 
ALPHAS(K)  =  TAUS(K)/TAUS(K-1) 

-  INTEGER  VARIABLE  USED  TO  CONTROL  ACCESS  TO  THE 
SUBROUTINE  ODD  WHEN  THE  INDEX  VARIABLE  K  IS  AN 
ODD  INTEGER. 

-  REAL  VECTOR  OF  "MODIFIED  NORM  VALUES".  THE 
VALUES  ARE  CALCULATED  FROM  A  SUMMATION  OF 
PRODUCT  TERMS  OF  THE  AUTOCORRELATION  LAGS,  AND 
THE  COEFFICIENTS  OF  THE  SINGULAR  PREDICTOR 
POLYNOMIALS. 

-  REAL  VECTOR  OF  REFLECTION  COEFFICIENTS  RHO(l), 
RHO(2),.  .  .  ,RHO(N). 

-  DESIRED  ORDER  OF  THE  PREDICTOR  POLYNOMIAL. 

-  REAL  VECTOR  OF  AUTOCORRELATION  LAGS  R(0),R(1), 
R(2),...  ,R(N) 

-  ARRAY  OF  SINGULAR  PREDICTOR  POLYNOMIAL 
COEFFICIENTS  FROM  P(0 ,0) , . . . P(N+1 ,N). 

-  ARRAY  OF  PREDICTOR  POLYNOMIAL  COEFFICIENTS. 
EACH  I-TH  ROW  OF  THE  ARRAY  CONTAINS  THE 
PREDICTOR  POLYNOMIAL  COEFFICIENTS  FOR  THE  I-TH 
ORDER  PREDICTOR  POLYNOMIAL. 

-  INTEGER  VARIABLE  USED  IN  THE  SUBROUTINE  ODD 
IN  THE  COMPUTATION  OF  THE  TAU(K)'S. 

-  DUMMY  VARIABLE  USED  DURING  THE  CALCULATION 
OF  THE  SYMMETRIC  SINGULAR  PREDICTOR  POLYNOMIAL 
COEFFICIENTS. 

VARIABLE  DECLARATIONS 


SIGMAN 
KEVEN 


LAMDA 
LAMDAS 


ALPHA 
ALPHAS 


KODD 


TAU 
TAUS 


RHO 
RHOS 
N 
C 

P 

PS 
A 
AS 
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REAL  C(0: 100) ,P(0: 100,0: 100) ,TAU(0: 100) 

REAL  A(0: 100,0: 100) , LAMDA( 0: 100) ,RHO( 100) , SIGMAN 

REAL  PS(0: 100,0: 100),AS(0: 100,0: 100) 

REAL  RHOS(100),ALPHAS(100),TAUS(0:  100),SIGMAS 

REAL  AR(0: 30) ,W(0: 5000) ,S(0:  5000) ,SIGMA(0:  100) 

REAL  R(0:  100) ,ALPHA( 100) , LAMDAS( 0:  100) 

REAL  AA(0: 100,0: 100) ,GAM( 50) ,LAMK,LAM( 100) 

INTEGER  M , LL , IX , T , KODD , KEVEN , L , N 

0PEN(L'NIT=4,BLANK='ZER0'  ) 
INITIALIZE  FILTER  ORDER 

READ(4,100)N 


LL  = 

N 

IX  = 

1913 

M  = 

3000 

WRITE(6; 

,300) 

WRITE(6: 

,400)N 

WRITE(6: 

,450)M 

DO  1 

1=0,  M 

CALL  GAUSS(IX,1.  ,0.  ,V) 

W(I)  =  V 
1       CONTINUE 

CALL  INPUT(LL,W,AR,S,M) 
C   INITIALIZE  AUTOCORRELATION  VECTOR  C 

CALL  ACORR(C,S,N+l,M) 

WRITE(6,2100) 
C    INITIAL  CONDITIONS  FOR  THE  SYMMETRIC  AND  ASYMMETRIC  PREDICTOR 
C   POLYNOMIAL  CALCULATIONS. 

P(0,0)  =  2. 

P(l,l)  =  I- 

TAU(O)  =  C(0) 
LAMDA(O)  =  1. 
KODD  =  1 
KEVEN  =  2 
A(N,0)  =  1.0 
PS(0,0)  =  0. 
PS(1,1)  =  -1. 
TAUS(O)  =  C(0) 
LAMDAS(O)  =  1. 
AS(N,0)  =  1.  0 

CALL  LEV(C, GAM, N,AA, LAM, SIGMA) 
WRITE(6,1700) 
WRITE(6,1800) 
DO  20  J=1,N 

WRITE(6,1900)N,J,SIGMA(J) ,LAM(J) ,AA(N,J) 
20     CONTINUE 
C   SYMMETRIC  &  ASYMMETRIC  SPLIT-LEVINSON  RECURSION 
WRITE(6,800) 
WRITE(6,850) 
WRITE (6, 9 00) 
DO  99  K=1,N 
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P(K,0)  =  1. 
PS(K,0)  =  1. 
C   CALL  STATEMENTS  FOR  EVEN  OR  ODD  VALUES  OF  INDEX  K 
IF(K.  EQ.  KODD)THEN 

CALL  AODD(C,PS,K,N,TAUS,T) 
CALL  ODD(C,P,K,N,TAU,KODD,T) 
ELSEIF(K.  EQ.  KEVEN)THEN 
CALL  AEVEN(C,PS,K,N,TAUS,T) 
CALL  EVEN(C,P,K,N,TAU,KEVEN,T) 
ENDIF 

ALPHA(K)  =  TAU(K)/TAU(K-1) 
ALPHAS(K)  =  TAUS(K)/TAUS(K-1) 
C   LOOP  TO  CALCULATE  SINGULAR  PREDICTOR  COEFFICIENTS 
DO  40  1=1, T 

P(K+1,I)  =  P(K,I)  +  P(K,I-1)  -  ALPHA(K)*P(K-1,I-1) 
PS(K+1,I)  =  PS(K,I)  +  PS(K,I-1)  -  ALPHAS(K)*PS(K-1,I-1) 
C   DECISION  PATH  TO  CALCULATE  SYMMETRIC  SINGULAR  PREDICTOR  COEFFICIENTS 
IFCI.LT.  T  .OR.  I.EQ.K)GOTO  40 
L  =  K+l 
DO  50  J=I+1,K 
P(L,J)  =  P(L,L-J) 
50       PS(L,J)  =  -PS(L,L-J) 
40     CONTINUE 

LAMDA(K)  =  2.  -  ( ALPHA(K)/LAMDA(K-1) ) 
LAMDAS(K)  =  2.  -  ( ALPHAS(K)/LAMDAS(K-1) ) 
C   REFLECTION  COEFFICIENT  CALCULATION 
RHO(K)  =  LAMDA(K)  -  1. 
RHOS(K)  =  1.  -  LAMDAS(K) 
WRITE(6, 1000)K,RHO(K) ,RHOS(K) ,GAM(K) 

99  CONTINUE 

C   CALCULATION  OF  N-TH  ORDER  NORM  (SIGMAN)  AND  N-TH  ORDER  PREDICTOR 
C   COEFFICIENTS,  A(N, 1) ,A(N,2) , . . . ,A(N,N) 

SIGMAN  =  LAMDA(N)*TAU(N) 

SIGMAS  =  LAMDAS(N)*TAUS(N) 

WRITE(6,1100) 

WRITE(6,600) 

WRITE(6,1200) 

DO  60  1=1, N 

A(N,I)  =  A(N,I-1)  +  P(N+1,I)  -  LAMDA(N)*P(N,I-1) 

AS(N,I)  =  AS(N,I-1)  +  PS(N+1,I)  -  LAMDAS(N)*PS(N,I-1) 

WRITE( 6 , 1300)1 ,TAU( I ) ,TAUS( I ) , ALPHA( I ) , ALPHAS( I ) , P( 1+1 , I ) , 
+PS(I+1,I),A(N,I),AS(N,I) 
60     CONTINUE 

100  F0RMAT(I2) 
200    F0RMAT(F12.6) 
300    FORMAT('l') 

400    FORMATC   FILTER  ORDER  =' ,13) 

450    FORMATC1-1,'  NUMBER  OF  SAMPLE  POINTS  =  ',15) 

600    FORMATC'-' ,103X,' FILTER  COEFFICIENTS') 

700    FORMAT(5X,I3,11X,F10.4,21X,F10.4) 
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800    FORMAT('-' ,21X, 'REFLECTION  COEFFICIENTS') 
850    FORMATC-'  ,25X,'SPLIT-LEVINSON') 

900    FORMAT(5X,' INDEX' , 8X, ' SYMMETRIC ' ,9X, ' ANTISYMMETRIC ' ,9X, 
+'LEVINSON' ) 
1000    FORMAT(5X,I3,10X,F12.  6,12X,F12.  6,12X,F12.  6) 

1100    FORMATC-'  ,5X, 'FILTER  PARAMETERS  FROM  SPLIT  LEVINSON  RECURSION') 
1200    FORMAT(5X,' INDEX' ,6X,'TAU(K)' ,4X, !TAU*(K) ' ,4X, ' ALPHA(K) ' ,4X, 

+'ALPHA*(K)' ,4X,'P(K)' ,4X,'P*(K)' ,6X, ' SYMMETRIC' ,6X, ' ASYMMETRIC' ) 
1300    FORMAT(5X,I3,6X,F12. 6,3X,F12. 6,3X,F12. 6,4X,F12.  6,3X,F12.  6,2X, 

+F12.  6,2X,F10.4,5X,F10.  4) 
1700    FORMATC-'  ,5X,  'FILTER  PARAMETERS  FROM  LEVINSON  RECURSION') 
1800   FORMATC-'  ,8x,'  INDEX'  ,  8X,  '  SIGMA(K)  '  ,5X,  '  LAMDA(K)  '  ,8X, 'FILTER 

+COEFFICIENTS' ) 
1900    FORMAT(8X,I2,I3,7X,F12. 6,6X,F12.  6,12X,F12.  6) 
2000    FORMAT(5X, ' SPLIT-LEVINSON  RECURSION  CALCULATIONS') 
2100    FORMAT(5X,' INDEX' ,5X, 'KNOWN  COEFFICIENTS' , 5X, ' AUTOCORRELATION 

+FUNCTION  C(K)') 
2200    F0RMAT(2X,I3,4X,F12. 6) 
2300    FORMAT(5X,I3,40X,F12.  6) 
WRITE(6,300) 
END 

SUBROUTINE  AODD(C,PS,K,N,TAUS ,T) 
C        THIS  SUBROUTINE  CALCULATES  THE  ANTISYMMETRIC  "MODIFIED 
C        NORMS"  WHEN  THE  INDEX  K  IS  AN  ODD  INTEGER. 
REAL  C(0: 100),PS(0: 100,0: 100),TAUS(0: 100) 
INTEGER  T 
T  =  (K-l)/2 
TEMP  =  0. 
DO  5  1=0, T 
5   TEMP  =  TEMP  +  (C(I)  -  C(K-I) )*PS(K, I) 
TAUS(K)  =  TEMP 
RETURN 
END 

SUBROUTINE  ODD( C , P , K , N , TAU , KODD , T) 
C        THIS  SUBROUTINE  CALCULATES  THE  "MODIFIED  NORMS"  (  TAU(K)'S) 
C        WHEN  THE  INDEX  K  IS  AN  ODD  INTEGER. 
REAL  C(0: 100), P(0: 100,0: 100),TAU(0: 100) 
INTEGER  T 
T  =  (K+l)/2 
TEMP  =  0. 
DO  15  1=0, T-l 
15   TEMP  =  TEMP  +  (C(I)  +  C(K-I))*P(K,I) 
TAU(K)  =  TEMP 
KODD  =  KODD  +  2 
RETURN 
END 

SUBROUTINE  AEVEN(C,PS ,K,N,TAUS ,T) 
C        SUBROUTINE  CALCULATES  THE  VALUE  OF  THE  ANTISYMMETRIC 
C        "MODIFIED  NORMS"  (TAUS(K)'S)  WHEN  THE  INDEX  K  IS  AN  EVEN 
C        INTEGER. 
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REAL  C(C  100), PS(0: 100,0: 100), TAUS(0: 100) 

INTEGER  r 

T  =  K/2 

TEMP=  0. 

PS(K,T)  =  0. 

DO  25  1=0, T-l 
25   TEMP  =  TEMP  +  (C(I)  -  C(K-I) )*PS(K, I) 

TAUSfK)  =  TEMP 

RETURN 

END 

SUBROUTINE  EVEN(C ,P,K,N,TAU,KEVEN,T) 
C        SUBROUTINE  CALCULATES  THE  VALUE  OF  THE  "MODIFIED  NORMS'* 
C        (TAU(K)'S)  WHEN  THE  INDEX  K  IS  AN  EVEN  INTEGER. 

REAL  C(0: 100), P(0: 100,0: 100),TAU(0: 100) 

INTEGER  T 

T  =  K/2 

TEMP=  0. 

DO  35  1=0, T-l 
35   TEMP  =  TEMP  +  (C(I)  +  C(K-I> )*P(K,I) 

TAU(K)  =  TEMP  +  C(T)*P(K,T; 

KEVEN  =  KEVEN  +  2 

RETURN 

END 

SUBROUTINE  INPUT(LL,W,AR,S ,M) 
C     SUBROUTINE  TO  GENERATE  THE  INPUT  SEQUENCE  FROM  A  GIVEN  FIR 
C     FILTER  AND  ZERO  MEAN,  UNIT  VARIANCE  WHITE  NOISE. 

REAL  W(0:M),AR(0:LL),S(0:M) 
C   CALCULATE  INPUT  SEQUENCE  VALUES  BETWEEN  0  AND  FILTER  ORDER. 
C     TEMP  =  0. 

S(0)  =  W(0) 

DO  45  K=1,M 

S(K)  =  W(K)  -0.6*W(K-1)  +  0.8*S(K-1) 
45   CONTINUE 

RETURN 

END 

SUBROUTINE  ACORR(C ,S ,NLAG,M) 
C     SUBROUTINE  TO  CALCULATE  THE  AUTO  CORRELATION  FUNCTION  OF  THE 
C     GIVEN  INPUT  SEQUENCE. 

REAL  C(0:  100)  ,S(0:  5000) 

INTEGER  T 

DO  5  I=0,NLAG 

TEMP  =  0. 

DO  15  T=0,M-1-I 

TEMP  =  TEMP  +  S(T)*S(T+I) 
15   CONTINUE 

C(I)  =  TEMP*(1.  /FLOAT  (M-l-I)) 
5   CONTINUE 

RETURN 

END 
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SUBROUTINE  LEV( C , GAM , N , AA , LAM , S IGMA ) 
C        SUBROUTINE  TO  CALCULATE  THE  PREDICTOR  FILTER  COEFFICIENTS 
C        USING  THE  LEVINSON  RECURSION. 

REAL  AA(0: 100,0: 100 ) ,C( 0: N+2) ,GAM( N) ,LAM( 100) 
REAL  LAMK,SIGMA(0:  100) 
C        INITIAL  CONDITIONS 
AA(0,0)  =  1. 
SIGMA(O)  =  C(0) 
C        COMPUTE  LAM(K),  RHO(K);  UPDATE  SIGMA(K)  FOR  THE  NEXT 
C        ITERATION. 
DO  10  K=1,N 
LAMK  =  0. 
AA(K,0)  =  1.  0 
DO  20  1=0, K-l 
LAMK  =  LAMK  -  C(K-I )*AA(K-1 , I) 
LAM(K)  =  LAMK 
20     CONTINUE 

GAM(K)  =  LAM(K)/SIGMA(K-1) 
SIGMA(K)  =  SIGMA(K-l)  -  LAM(K)*GAM(K) 
C        COMPUTE  THE  PREDICTOR  FILTER  COEFFICIENTS 
IF(K  .EQ.  1)THEN 
AA(1,1)  =  GAM(K) 
ELSE 
DO  30  1=1, K-l 
AA(K,I)  =  AA(K-1,I)  +  GAM(K)*AA(K-1,K-I) 
30      CONTINUE 

AA(K,K)  =  GAM(K) 
ENDIF 
10     CONTINUE 
RETURN 
END 


4') 


APPENDIX  C.      SPLIT  LATTICE  ALGORITHMS 

C  PROGRAM  TO  CALCULATE  THE  NTH  ORDER  LATTICE  REFLECTION 

C  COEFFICIENTS  FROM  A  GIVEN  SEQUENCE  USING  THE  SYMMETRIC  ERROR 

C  VECTOR,  THE  ASYMMTRIC  ERROR  VECTOR,  OR  THE  FORWARD  AND 

C  BACKWARD  ERROR  VECTORS.  VARIABLES  DEFINED  IN  PREVIOUS  APPENDICES 

C  ARE  NOT  REDEFINED. 

C  *****  THIS  PROGRAM  REQUIRES  SUBROUTINE  INPUT  FROM  APPENDIX  A  ******* 

C  VARIABLE  DEFINITIONS 

c  SIG    -  N-TH  DEGREE  NORM  OF  THE  FILTER. 

C  GAM    -  VECTOR  OF  REFLECTION  COEFFICIENTS  CALCULATED  BY 

C  THE  LEVINSON  RECURSION. 

C  LAM    -  REAL  VARIABLE  USED  WHEN  CALCULATING  THE  REFLECTION 

C  COEFFICIENTS  FROM  THE  LEVINSON  RECURSION 

C  THE  REFLECTION  COEFFICIENT  IN  TERMS 

C  OF  THE  FILTER  NORM  IS  GIVEN  BY: 

C  RHO(K)  =  LAM/SIG 

C  TAU    -  REAL  VECTOR  OF  "MODIFIED  NORM  VALUES".  THE 

C  VALUES  ARE  CALCULATED  FROM  A  SUMMATION  OF 

C  PRODUCT  TERMS  OF  THE  SYMMETRIC  OR  ASYMMETRIC 

C  PREDICTION  ERROR  SEQUENCES. 

C  A      -  ARRAY  OF  PREDICTOR  POLYNOMIAL  COEFFICIENTS. 

C  AS       EACH  I-TH  ROW  OF  THE  ARRAY  CONTAINS  THE 

C  AL       PREDICTOR  POLYNOMIAL  COEFFICIENTS  FOR  THE  I-TH 

C  ORDER  PREDICTOR  POLYNOMIAL. 

C  AR       VECTOR  OF  COEFFICIENTS  FROM  THE  KNOWN  TEST  FILTER. 

C  XO     -  SYMMETRIC  OR  ASSYMMETRIC  PREDICTION  ERROR  VECTOR 

C  FOR  THE  (K-l)  STAGE  OF  THE  LATTICE  FILTER. 

C  LL     -  DESIRED  LATTICE  FILTER  ORDER. 

C  XI      -  SYMMETRIC  OR  ASYMMETRIC  PREDICTION  ERROR  VECTOR 

C  FOR  THE  K-TH  STAGE  OF  THE  LATTICE  FILTER. 

C  AT     -  TEMP  STORAGE  FOR  THE  PREDICTION  ERROR  VECTOR  WHILE 

C  COMPUTING  THE  (K+l)  STAGE  PREDICTION  ERROR  VECTOR. 

C  FT     -  SHIFTED  FORWARD  PREDICTION  ERROR  VECTOR. 

C  BT     -  SHIFTED  BACKWARD  PREDICTION  EROR  VECTOR. 

C  M      -  DESIRED  ORDER  OF  THE  PREDICTOR  POLYNOMIAL. 

C  T      -  INTEGER  VARIABLE  USED  IN  THE  PROGRAM. 

C  W      -  WHITE  NOISE  SEQUENCE  VECTOR. 

C  S      -  INPUT  SEQUENCE  VECTOR 

C  F        FORWARD  PREDICTION  ERROR  VECTOR. 

C  B        BACKWARD  PREDICTION  ERROR  VECTOR. 

C  VARIABLE  DECLARATIONS 

REAL  AR(30) ,W(0: 5000) ,S(0: 5000) ,RHO(100) 

REAL  A(0:  100,0:  100) ,GAM(20) ,RHOS(  100) ,AS(0:  100,0:  100) 

REAL  ALPHA, X1(0: 5000), XO(0:  5000), AT(0:  5000), AL(0:  100,0: 100) 

INTEGER  M,LL,IX,T,L,N 

0PEN(UNIT=4,BLANK='ZER0*) 

C  INITIALIZE  FILTER  ORDER 
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READ(4,100)M 
C    INITIAL  CONDITIONS  FOR  INPUT  SEQUENCE  GENERATOR 
LL  =  M 
N  =  5000 
IX  =  1913 
WRITE(6,200) 
WRITE(6,300)M 
WRITE(6,400)N 
WRITE(6,500) 
WRITE(6,600) 
C   LOOP  TO  GENERATE  WHITE  NOISE  SEQUENCE  AND  TO  READ  TEST  COEFFICIENTS. 
DO  1  1=0, N 
CALL  GAUSS(IX,1.  ,0.  ,V) 
W(I)  =  V 
1      CONTINUE 
C   CALL  STATEMENT  TO  GENERATE  INPUT  SEQUENCE 

CALL  INPUT(LL,W,AR,S,N) 
C        CALL  STATEMENTS  FOR  SYMMETRIC,  ASYMMETRIC,  AND  LEVINSON  LATTICE 
C        RECURSION  SUBROUTINES 

CALL  SLAT(S,M,N,RHO,ALPHA,Xl,AT,XO) 
CALL  ASLAT(S,M,N,RHOS,ALPHA,Xl,AT,XO) 
CALL  LEV(N,GAM,M,S) 
WRITE(6,700) 
WRITE(6,800) 
DO  80  K=1,M 

WRITE(6,900)K,RHO(K),RHOS(K),GAM(K) 
80     CONTINUE 

DO  90  K=1,M 
A(K,0)  =  1. 
AS(K,0)  =  1. 
AL(K,0)  =  1. 
IF(K  .EQ.  1)THEN 
A(l,l)  =  RHO(K) 
AS(1,1)  =  RHOS(K) 
AL(1,1)  =  GAM(K) 
ELSE 
DO  95  1=1, K-l 
A(K,I)  =  A(K-1,I)  +  RHO(K)*A(K-l,K-I) 
AS(K,I)  =  AS(K-1,I)  +  RHOS(K)*AS(K-l,K-I) 
AL(K,I)  =  AL(K-1,I)  +  GAM(K)*AL(K-1,K-I) 
95       CONTINUE 

A(K,K)  =  RHO(K) 
AS(K,K)  =  RHOS(K) 
AL(K,K)  =  GAM(K) 
ENDIF 
90     CONTINUE 

WRITE(6,1000) 
WRITE(6,1100) 
WRITE(6,1200) 
DO  96  K=1,M 


WRITE(6,1300)M,K,AL(M,K),A(M,K),AS(M,K) 
96     CONTINUE 
100    F0RMAT(I2) 
200    FORMAT('l') 

300    FORMATC   FILTER  ORDER  =  ',13) 
400    FORMATC  V   NUMBER  OF  SAMPLE  POINTS  =  ',15) 
500    FORMATC '-' ,10X, 'KNOWN  FILTER  COEFFICIENTS' ) 
600    FORMATC-' ,8X,' INDEX' ,10X, 'FILTER  COEFFICIENTS') 
700    FORMATC-' ,10X,' REFLECTION  COEFFICIENTS*) 
800    FORMATC'-' ,5X,'  INDEX  ',  3X,  '  SYMMETRIC'  ,  9X,  *  ANTISYMMETRIC  ' 

+,9X,'  LEVINSON  ') 
900    FORMATC-'  ,6X,  13  ,6X,F8.  4, 12X.F8.  4, 12X.F8.  4) 
1000    FORMATC'-' ,15X,'  FILTER  COEFFICIENTS  ') 
1100    FORMATC'-' ,20X,'  LEVINSON  ' ,12X,'  SPLIT-LEVINSON  *) 
1200    FORMATC 5X,'  INDEX  ',26X,'  SYMMETRIC  ',4X,'  ASYMMETRIC  ') 
1300    FORMATC'  ' , 6X, 212 , 10X,F8.  4, 11X,F8. 4, 7X,F8. 4) 
WRITE(6,200) 
END 

SUBROUTINE  SLAT(S,M,N,RHO, ALPHA, XI , AT,XO) 
C        SUBROUTINE  TO  COMPUTE  THE  LATTICE  REFLECTION  COEFFICIENTS 
C        USING  THE  SYMMETRIC  PREDICTION  ERROR  VECTOR. 
REAL  X1(0: M+N) ,X0C0: M+N) ,RHO(M) ,SC0:  N) , ALPHA 
REAL  AT(0:M+N),A(100,100) 
C    INITIAL  CONDITIONS 
INTEGER  T 
RRHO  =  0. 
TAU  =  0. 
C   INITIALIZE  THE  PREDICTION  ERROR  VECTORS  FOR  THE  ZERO  AND  1ST 
C   STAGES  OF  THE  LATTICE.  INITIALIZE  THE  ZERO  STAGE  MODIFIED  NORM 
DO  10  T=0,N-1 
XO(T)  =  2.*S(T) 
TAU  =  TAU  +  S(T)**2 
10     CONTINUE 

DO  20  T=0,N 
IF(T.  EQ.  N)S(T)  =  0. 
IF(T.  EQ.  0)THEN 

X1(T)  =  SCT) 
ELSE 

X1(T)  =  S(T)  +  S(T-l) 
ENDIF 
20     CONTINUE 
C    LOOP  TO  COMPUTE  THE  REFLECTION  COEFFICIENTS 
DO  101  K=1,M 
TTAU  =  TAU 
C     STORE  TAU(K-l),  AND  COMPUTE  TAU(K). 
TAU  =  0. 
DO  30  T=0,N+K-1 
TAU  =  TAU  +  X1(T)**2 
30     CONTINUE 
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TAU  =  TAU/2. 
C     COMPUTE  ALPHA(K),  RHO(K);  STORE  TAU(K)  AND  RHO(K)  FOR 
C     NEXT  ITERATION. 

ALPHA  =  TAU/TTAU 
TTAU  =  TAU 

RHO(K)  =  1.  -  ( ALPHA/ (1.  +  RRHO)) 
RRHO  =  RHO(K) 
C     LOOP  TO  COMPUTE  THE  (K+l)  PREDICTION  ERROR  VECTOR. 
DO  40  T=0,N+K 
IF(T  .  EQ.  0)THEN 
AT(T)  =  X1(T) 
ELSEIF(T.  EQ.N+K)THEN 
AT(T)  =  Xl(T-l) 
ELSE 

AT(T)  =  X1(T)  +  Xl(T-l)  -  ALPHA*XO(T-l) 
ENDIF 
40     CONTINUE 
C     LOOPS  TO  UPDATE  PREDICTION  ERROR  VECTORS  FOR  NEXT  ITERATION. 
C         (SHIFT  XI  INTO  XO  AND  AT  INTO  XI) 
DO  50  T=0,N+K-1 
XO(T)  =  X1(T) 
50     CONTINUE 

DO  60  T=0,N+K 
X1(T)  =  AT(T) 
60     CONTINUE 
101    CONTINUE 
RETURN 
END 

SUBROUTINE  ASLAT(S ,M,N,RHOS .ALPHA, XI , AT.XO) 
C        SUBROUTINE  TO  COMPUTE  THE  LATTICE  REFLECTION  COEFFICIENTS 
C        USING  THE  ASYMMETRIC  PREDICTION  ERROR  VECTOR. 
REAL  Xl( 0: M+N) ,XO(0: M+N) ,RHOS(M) ,S(0:  N) , ALPHA 
REAL  AT(0:M+N) 
INTEGER  T 
C    INITIAL  CONDITIONS 
RRHO  =  0. 
TAU  =  0. 
C   INITIALIZE  THE  PREDICTION  ERROR  VECTORS  FOR  THE  ZERO  AND  1ST 
C   STAGES  OF  THE  LATTICE.  INITIALIZE  THE  ZERO  STAGE  MODIFIED  NORM 
DO  10  T=0,N-1 
XO(T)  =  0. 

TAU  =  TAU  +  S(T)**2 
10     CONTINUE 

DO  20  T=0,N 
IF(T.  EQ.  N)X1(T)  =  -S(T-l) 
IF(T.  EQ.  0)THEN 

X1(T)  =  S(T) 
ELSE 

X1(T)  =  S(T)  -  S(T-l) 
ENDIF 
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20     CONTINUE 
C    LOOP  TO  COMPUTE  THE  REFLECTION  COEFFICIENTS 

DO  101  K=1,M 
C     STORE  TAU(K-l),  AND  COMPUTE  TAU(K). 
TTAU  =  TAU 
TAU  =  0. 
DO  30  T=0,N+K-1 
TAU  =  TAU  +  X1(T)**2 
30     CONTINUE 

TAU  =  TAU/ 2. 
C     COMPUTE  ALPHA(K),  RHO(K);  STORE  TAU(K)  AND  RHO(K)  FOR 
C     NEXT  ITERATION. 

ALPHA  =  TAU/TTAU 
TTAU  =  TAU 

RHOS(K)  =  ( ALPHA/ (1.  -  RRHO))  -  1. 
RRHO  =  RHOS(K) 
C     LOOP  TO  COMPUTE  THE  (K+l)  PREDICTION  ERROR  VECTOR. 
DO  40  T=0,N+K 
IF(T  .EQ.  0)THEN 

AT(T)  =  X1(T) 
ELSEIF(T. EQ.  N+K)THEN 
AT(T)  =  Xl(T-l) 
ELSE 

AT(T)  =  X1(T)  +  Xl(T-l)  -  ALPHA*X0(T-1) 
ENDIF 
40     CONTINUE 
C     LOOPS  TO  UPDATE  PREDICTION  ERROR  VECTORS  FOR  NEXT  ITERATION. 
C         (SHIFT  XI  INTO  XO  AND  AT  INTO  XI) 
DO  50  T=0,N+K-1 
XO(T)  =  X1(T) 
50     CONTINUE 

DO  60  T=0,N+K 
X1(T)  =  AT(T) 
60     CONTINUE 
101    CONTINUE 
RETURN 
END 

SUBROUTINE  LEV(N,GAM,M,S) 
C        SUBROUTINE  TO  COMPUTE  THE  REFLECTION  COEFFICIENTS  FOR  AN 
C        N-TH  ORDER  LATTICE  FILTER  FROM  THE  FORWARD  AND  BACKWARD 
C        PREDICTION  ERROR  VECTORS. 

REAL  F(0:  5100), B(0:  5100), FT(0: 5100), BT(0: 5100), GAM(20) 
REAL  LAM,SIG,S(0:N) 
INTEGER  T 
C     INITIAL  CONDITIONS;  INITIALIZE  FORWARD  AND  BACKWARD  PREDICTION 
C     ERROR  VECTORS. 
SIG  =  0. 
DO  10  1=0, N-l 
F(I)  =S(I) 
B(I)  =  S(I) 
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ftc  i)  =  sen 

SIG  =  SIG  +  S(I)**2 

10     CONTINUE 
C     LOOP  TO  COMPUTE  THE  REFLECTION  COEFFICIENTS 

DO  20  K=1,M 
C     FORM  SHIFTED  ERROR  VECTORS 
DO  30  T=1,N+K-1 
BT(T)  =  B(T-l) 
30     CONTINUE 
BT(0)  =  0. 
FT(N+K-1)  =  0. 
LAM  =  0. 
C     COMPUTE  LAM(K),  GAM(K);  UPDATE  K-TH  ERROR  NORM  AND 
C     STORE  FOR  NEXT  ITERATION. 
DO  40  T=l,N+K-2 
LAM  =  LAM  -  FT(T)*BT(T) 
40     CONTINUE 

GAM(K)  =  LAM/ SIG 
IF(K  .EQ.  M)GOTO  20 
SIG  =  SIG  -  LAM*GAM(K) 
C     COMPUTE  (K+l)  FORWARD  AND  BACKWARD  PREDICTION  ERRORS  AND  SHIFT 
C     INTO  TEMPORARY  VECTORS  FOR  NEXT  ITERATION. 
DO  50  T=0,N+K-1 
F(T)  =  FT(T)  +  GAM(K)*BT(T) 
B(T)  =  GAM(K)*FT(T)  +  BT(T) 
50     CONTINUE 

DO  60  T=0,N+K-1 
FT(T)  =  F(T) 
BT(T)  =  B(T) 
60     CONTINUE 
20     CONTINUE 
RETURN 
END 
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APPENDIX  D.      MA  PREDICTOR  COEFFICIENT  PROGRAM 


THIS  PROGRAM  IS  TO  COMPUTE  THE  FIR  FILTER  COEFFICIENTS 
USING  THE  SPLIT-LEVINSON  ALGORITHM,  AND  AN  AUTOREGRESSIVE 
MOVING  AVERAGE  PROCESS.  VARIABLES  DEFINED  IN  PREVIOUS 
APPENDICES  ARE  NOT  REDEFINED. 
THIS  PROGRAM  REQUIRES  THE  SUBROUTINES  ODD,  AND  EVEN 
FROM  APPENDIX  A  **** 

VARIABLE  DEFINITIONS 

-  PREDICTOR  COEFFICIENT  ERROR  NORM. 
ENORM  =  SQRT((A(0)-AA(0))**2  +. . . +  (A(N)-AA(N) )**2) 

-  NUMBER  OF  POINTS  OF  INPUT  SEQUENCE  TO  START. 

-  NUMBER  OF  POINTS  OF  INPUT  SEQUENCE  AT  END  OF  PROGRAM. 

-  ERROR  VECTOR. 

-  ERROR  VECTOR. 

-  NUMBER  OF  INPUT  DATA  POINTS  (0,1,... ,NPTS). 

-  BACKWARD  PREDICTION  ERROR  POWER  OF  X. 

-  VECTOR  OF  X  AND  Y  CROSSCORRLATION  ELEMENTS 

-  VECTOR  OF  CALCULATED  PREDICTOR  COEFFICIENTS. 

-  INDEX  FOR  X-AXIS  OF  PLOTTING  ROUTINE. 

-  VECTOR  OF  PREDICTOR  COEFFICIENT  NORMS. 

-  FORWARD  PREDICTION  ERROR  POWER  OF  X. 

-  FORWARD  PREDICTION  ERROR  POWER  OF  Y. 

-  Y  REFLECTION  COEFFICIENT. 

-  X  REFLECTION  COEFFICIENT. 

-  FILTER  ORDER  VARIABLE  USED  IN  SUBROUTINE  CORR. 

-  INTEGER  SEED  NUMBER  USED  BY  IMSL  SUBROUTINE  GAUSS. 

-  X  AUTOCORRELATION  VECTOR. 

-  Y  AUTOCORRELATION  VECTOR. 

-  MA  COEFFICIENT  VECTOR. 

-  MA  COEFFICIENT  VECTOR. 

-  MA  COEFFICIENT  VECTOR. 

-  INTEGER  VARIABLE  USED  IN  THE  SUBROUTINE  ODD 
IN  THE  COMPUTATION  OF  THE  TAU(K)'S. 

L      -  DUMMY  VARIABLE  USED  DURING  THE  CALCULATION 

OF  THE  SYMMETRIC  SINGULAR  PREDICTOR  POLYNOMIAL 
COEFFICIENTS. 
X      -  INPUT  WHITE  NOISE  VECTOR. 
M      -  INDEXING  VARIABLE  USED  IN  FIR  FILTER  COEFFICIENT 

RECURSION  (M=l,2,. . . ,N). 
Y      -  OUTPUT  SEQUENCE  FROM  FIR  TEST  FILTER. 
VARIABLE  DECLARATIONS 
REAL  P(0: 100,0: 100),TAU(0: 100) ,C(0: 50) ,B(0: 50) ,EN(200) 
REAL  A(0:  100,0:  100) ,LAMDA( 0:  100), X(0:  10000)  ,D(0: 50) 
REAL  AR(0: 30),Y(0: 10000) ,EY( 0: 50) ,EX( 0: 50) ,KX( 50) 
REAL  DELX(0: 50) ,DELY(0: 50) ,EXB(0: 50) ,KY(50) ,XX(200) 
REAL  RX(0:  100) ,ALPHA( 100) ,RY(0: 2) ,RXY(0:  100),AA(0:  100) 
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INTEGER  M , LL , I X , T , KODD , KE  VEN , L , N , NPTS 
C    DESIRED  FILTER  ORDER  AND  THE  TEST  FILTER  COEFFICIENTS 
C         ARE  READ  FROM  A  DATA  FILE  (FILE  FT04F001). 

OPEN( UNIT=4 , BLANK= ' ZERO ' ) 
C   INITIALIZE  FILTER  ORDER  AND  TEST  FILTER  COEFFICIENTS 

READ(4,100)N 

LL  =  N 

DO  20  K=0,LL 

READ(4,700)AR(K) 
20      CONTINUE 

NPTS  =  4000 
WRITE(6,300) 
WRITE(6,400)N 

IX  =  1913 
WRITE(6,600)NPTS 
C    GENERATE  (NPTS+1)  WHITE  NOISE  SEQUENCE 

DO  10  1=0, NPTS 

CALL  GAUSS(IX,1.  ,0.  ,V) 

X(I)  =  V 
10      CONTINUE 
C   CREATE  INPUT  SEQUENCE  FROM  GIVEN  FIR  TEST  FILTER 

CALL  INPUT(LL,X,AR,Y,NPTS) 
C   INITIALIZE  AUTOCORRELATION  VECTOR  RX,RY,  AND  CROSSCORRELATION 
C   VECTOR  RXY 

CALL  CORR( N+l, NPTS, X,Y,RX,RY, RXY) 

WRITE(6,800) 

DO  30  1=0, N 

WRITE( 6 , 900) I , AR( I ) ,RX( I ) 
30     CONTINUE 
C   INITIAL  CONDITIONS  FOR  MOVING  AVERAGE  MODEL  PARAMETERS 

BOO  =  -RXY(0)/RX(0) 

EY(0)  =  RY(0)  -  B00*RXY(0) 

EX(0)  =  RX(0) 

DO  40  1=0,1 

C(I)  =  I 

D(I)  =  I 
40     CONTINUE 

B(0)  =  1. 

B(l)  =  BOO 

DELY(O)  =  RXY(l)  -  B00*RX(1) 

DELX(O)  =  RX(1) 
C   LOOP  TO  GENERATE  PREDICTOR  COEFFICIENT  VECTOR 

DO  120  M=1,N 
C   INITIAL  CONDITIONS  FOR  THE  SYMMETRIC  PREDICTOR  POLYNOMIAL 
C   CALCULATIONS. 

P(0,0)  =  2. 

P(1,D  =  1. 

TAU(O)  =  RX(O) 

LAMDA(O)  =  1. 

KODD  =  1 
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KEVEN  =  2 
A(M,0)  =  1. 
C   SYMMETRIC  SPLIT-LEVINSON  RECURSION 
DO  70  K=1,M 
P(K,0)  =  1. 
C   CALL  STATEMENTS  FOR  EVEN  OR  ODD  VALUES  OF  INDEX  K 
IF(K.  EQ.  KODD)THEN 

CALL  ODD(RX,P,K,N,TAU,KODD,T) 
ELSEIF(K.  EQ.  KEVEN)THEN 
CALL  EVEN(RX,P,K,N,TAU,KEVEN,T) 
END  IF 

ALPHA(K)  =  TAU(K)/TAU(K-1) 
C   LOOP  TO  CALCULATE  SINGULAR  PREDICTOR  COEFFICIENTS 
DO  60  1=1, T 

P(K+1,I)  =  P(K,I)  +  P(K,I-1)  -  ALPHA(K)*P(K-1,I-1) 
C   DECISION  PATH  TO  CALCULATE  SYMMETRIC  SINGULAR  PREDICTOR  COEFFICIENTS 
IF(I.LT.  T  .OR.  I.EQ.  K)GOTO  60 
L  =  K+l 
DO  50  J=I+1,K 
P(L,J)  =  P(L,L-J) 
50      CONTINUE 
60     CONTINUE 

LAMDA(K)  =  2.  -  ( ALPHA(K)/LAMDA(K-1) ) 
70     CONTINUE 
C   CALCULATION  OF  N-TH  ORDER  PREDICTOR  COEFFICIENTS,  A(K,  1) ,.  .  .  A(K,K) 
C      COMPUTE  PREDICTOR  COEFFICIENTS  FOR  K-TH  ITERATION 
DO  80  1=1, M 
A(M,I)  =  A(M,I-1)  +  P(M+1,I)  -  LAMDA(M)*P(M,I-1) 
C(I)  =  A(M,I) 
80     CONTINUE 

DO  90  J=1,M 
D(l)  =  -C(J) 
IF(J  .LT.  M)THEN 
DO  95  I=J+1,M 
D(I)  =  D(I-l) 
95       CONTINUE 

ENDIF 
90     CONTINUE 
D(0)  =  0. 
D(M+1)  =  1. 
C   UPDATE  PREDICTION  ERRORS 
XBTMP  =  0. 
XTMP  =  0. 
DO  25  1=1, M 
XBTMP  =  XBTMP  +  C( I)*RXY(M+1-I) 
XTMP  =  XTMP  +  C(I)*RX(M+1-I) 
25     CONTINUE 

DELX(M)  =  RX(M+1)  -  XTMP 
EXB(M)  =  RXY(M+1)  -  XBTMP 
C   UPDATE  REFLECTION  COEFFICIENTS 
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KX(M)  =  -DELX(M-1)/EX(M-1) 
EX(M)  =  EX(M-l)  +  KX(M)*DELX(M-1) 
KY(M)  =  -DELY(M-1)/EX(M) 
EY(M)  =  EY(M-l)  +  KY(M)*EXB(M) 
C   UPDATE  B  VECTOR 
B(M+1)  =  0. 
DO  45  1=0, M+l 
B(I)  =  B(I)  +  KY(M)*D(I) 
45     CONTINUE 
YTMP  =  0. 
DO  55  1=1, M+l 

YTMP  =  YTMP  -  B(I)*RX(M+2-I) 
55     CONTINUE 

DELY(M)  =  RXY(M+1)  -  YTMP 
IF(M  . EQ.  N)THEN 
WRITE(6,1000)N 
WRITE(6,1100) 
DO  130  K=0,M 
WRITE(6,1200)K,-B(K+1) 
AA(K)  =  -B(K+1) 
130     CONTINUE 

ENDIF 
120    CONTINUE 
100    FORMAT(I2) 
200    FORMAT( '  ' ) 
300    FORMAT('l') 

400    FORMATC   FILTER  ORDER  =' ,13) 
500    FORMATC1-') 

600    FORMATC -V  NUMBER  OF  SAMPLE  POINTS  =  ',15) 
700    FORMATC F8.  4) 
800    FORMATC'-' ,'5X,' INDEX' ,5X,' KNOWN  COEFFICIENTS  * ,5X, 

+' AUTOCORRELATION  FUNCTION  R(K)') 
900    FORMATC  '  ,5X, 13 , 11X.F8.  4 , 21X,F8.  4) 
1000    FORMATC'-' ,10X,' PREDICTOR  COEFFICIENTS  FOR  FILTER  ORDER  =  ',13) 
1100    FORMATC'-' ,5X,' INDEX' ,12X,' FIR  PREDICTOR  COEFFICIENTS') 
1200    FORMATC  '  ,5X,  13  ,  23X,F8.  4) 
WRITE(6,300) 

END 
SUBROUTINE  INPUT(LL,X,AR,Y,NPTS) 
C     SUBROUTINE  TO  GENERATE  THE  INPUT  SEQUENCE  FROM  A  GIVEN  FIR 
C     FILTER  AND  ZERO  MEAN,  UNIT  VARIANCE  WHITE  NOISE. 

REAL  X(0:NPTS),AR(0:LL),Y(0:NPTS) 
C   CALCULATE  INPUT  SEQUENCE  VALUES  BETWEEN  0  AND  FILTER  ORDER. 
DO  45  K=0,NPTS 
TEMP  =  0. 
IF(K.  LE.LDTHEN 

LU=K 
ELSE 

LU=LL 
ENDIF 
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J  =  K 

DO  55  1=0, LU 
TEMP  =  TEMP  +  AR(I)*X(J) 
J  =  J-l 
55   CONTINUE 

Y(K)  =  TEMP 
45   CONTINUE 
RETURN 
END 

SUBROUTINE  CORR(NLAG,NPTS ,X, Y,RX,RY,RXY) 
C     SUBROUTINE  TO  CALCULATE  THE  AUTOCORRELATION  FUNCTION  X  AND  Y  AND 
C     THE  CROSSCORRELATION  FUNCTION  BETWEEN  X  AND  Y 

REAL  RX(0: NLAG) ,Y(0: NPTS) ,X(0:  NPTS) ,RXY(0:  NLAG) ,RY(0:  2) 
INTEGER  T 
C   COMPUTE  THE  AUTOCORRELATION  OF  X  AND  THE  CROSSCORRELATION  OF 
C   X  AND  Y  FOR  LAGS  0  , 1 ,  2 , . . .  ,NLAG 
DO  5  1=0, NLAG 
XTEMP  =  0. 
XYTEMP  =  0. 
C   COMPUTE  THE  AUTOCORRELATION  OF  X  AND  THE  CROSSCORRELATION  OF 
C   X  &  Y  FOR  LAG  I 

DO  15  T=0,NPTS-1-I 
XTEMP  =  XTEMP  +  X(T)*X(T+I) 
XYTEMP  =  XYTEMP  +  X(T)*Y(T+I) 

15  CONTINUE 

RX(I)  =  XTEMP*(1.  /FLOAT(NPTS-l-I)) 
RXY(I)  =  XYTEMP*(1./FL0AT(NPTS-1-I)) 
C   COMPUTE  THE  ZERO  LAG  AUTOCORRELATION  FUNCTION  OF  Y 
IF(I  .EQ.  0)THEN 
RY(O)  =  0. 
DO  16  J=0,NPTS-1 
RY(O)  =  RY(O)  +  Y(J)**2 

16  CONTINUE 

RY(O)  =  RY(0)*(1./FLOAT(NPTS-1)) 
END  IF 
5      CONTINUE 
RETURN 
END 
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APPENDIX  E.     EXTENDED  PRONY  PROGRAM 

C        PROGRAM  TO  CALCULATE  THE  NTH  ORDER  LATTICE  REFLECTION 

C        COEFFICIENTS  FROM  A  GIVEN  SEQUENCE  USING  THE  SYMMETRIC  ERROR 

C        VECTOR,  THE  ASYMMTRIC  ERROR  VECTOR,  OR  THE  FORWARD  AND 

C        BACKWARD  ERROR  VECTORS. 

C  VARIABLE  DEFINITIONS 

C         PT   -  TEMPORARY  ARRAY  USED  TO  AVERAGE  PREDICTOR 

C  COEFFICIENTS. 

C         PP    -  ESTIMATED  NUMBER  OF  COMPLEX  SINUSOIDS  PRESENT. 

C         Al    -  AMPLITUDE  OF  #1  SINUSOID,  (1-4)  SINUSOIDS 

C  PRESENT. 

C         FS    -  SAMPLING  FREQUENCY. 

C         Fl    -  FREQUENCY  OF  #1  SINUSOID  IN  TEST  SEQUENCE. 

C      THETA1   -  DIGITAL  FREQUENCY  OF  #1  TEST  ANALOG  FREQUENCY. 

C  VARIABLE  DECLARATIONS 

REAL  W(0: 5000) ,S(0: 5000) ,ALPHA( 100) ,ROOTR( 100) ,XCOF(0:  100) 
REAL  P(0:  100,0:  100) , ALPHAS( 100) ,COF( 0:  100) ,ROOTI( 100) 
REAL  X1(0:5000),XO(0:  5000), AT(0:  5000), PS(0:  100,0: 100),PT(0: 100) 
INTEGER  T,PP 

OPEN( UNIT=4 , BLANK=' ZERO ' ) 
C   INITIALIZE  FILTER  ORDER 

READ(4,100)PP 

M  =  2*PP 
C   INITIAL  CONDITIONS  FOR  INPUT  SEQUENCE  GENERATOR 

IX  =  1913 

A  =  SQRT(2. ) 

B  =  SQRT(2.  ) 

C  =  SQRT( 10.  ) 
C  D  =  SQRT(2.0) 
C        E  =  SQRT(2.0) 

Fl  =  5.5E+01 

F2  =  7.5E+02 
C        F3  =  1.  25E+02 
C        F4  =  1.  75E+02 

FS  =  2.25E+02 

TWOPI  =  6.2831853 

THETA1  =  (TW0PI*F1)/FS 

THETA2  =  (TWOPI*F2)/FS 
C   LOOP  TO  GENERATE  WHITE  NOISE  SEQUENCE  AND  TO  READ  TEST  COEFFICIENT 

DO  1  1=0, N 

CALL  GAUSS(IX,1.  ,0.  ,V) 

W(I)  =  C*V 

Al  =  A*C0S(TW0PI*(F1/FS)*FL0AT(I)) 

A2  =  B*COS(TWOPI*(F2/FS)*FLOAT(I)) 
C  A3  =  D*COS(TWOPI*(F3/FS)*FLOAT(I)) 
C        A4  =  E*C0S(TW0PI*(F4/FS)*FL0AT(I)) 

S(I)  =  Al  +  A2  +  W(I) 
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I  CONTINUE 

C  CALL  STATEMENTS  FOR  SYMMETRIC,  ASYMMETRIC,  AND  LEVINSON  LATTICE 
C   RECURSION  SUBROUTINES 

CALL  SLAT(S,M,N,P,ALPHA,Xl,AT,XO) 
C        CALL  ASLAT(S,M,N,PS,ALPHAS,Xl,AT,XO) 
WRITE(6,200) 

WRITE(6,300)PP 
WRITE(6,400)M 
WRITE(6,500)N 
C        DISPLAY  COEFFICIENTS  OF  POLYNOMIAL 
WRITE(6,600) 
DO  11  K=0,M 
IF(K  .EQ.  M)P(M,K)  =  1.0 
WRITE(6,700)M,K,P(M,K) 

II  CONTINUE 
100    F0RMAT(I4) 
200    FORMAT('l') 

300    FORMATC  NUMBER  OF  COMPLEX  EXPONENTIALS  IN  MODEL  =  ',13) 
400    FORMATC '  V   SYMMETRIC  FILTER  ORDER  =  ',13) 
500    FORMATC  ','   NUMBER  OF  SAMPLE  POINTS  =  ',15) 
600    FORMATC-'  ,8X,' INDEX'  ,  13X, '  COEFFICIENTS'  ) 
700    FORMAT(5X,2I2,16X,F8.4) 
WRITE(6,200) 
END 
SUBROUTINE  SLAT(S,M,N,P,ALPHA,X1,AT,X0) 
C        SUBROUTINE  TO  COMPUTE  THE  SYMMETRIC  PREDICTOR  COEFFICIENTS 
C        USING  THE  SYMMETRIC  PREDICTION  ERROR  VECTOR. 
REAL  X1(0:  M+N) ,XO(0:  M+N) ,ALPHA(M) ,S(0: N) 
REAL  AT(0:  M+N), P(0:  100,0:  100) 
INTEGER  T 
C        INITIAL  CONDITIONS 
KODD  =  1 
KEVEN  =  2 
TAU  =  0. 
P(l.l)  =  1. 
P(0,0)  =  2. 
C   INITIALIZE  THE  PREDICTION  ERROR  VECTORS  FOR  THE  ZERO  AND  1ST 
C   STAGES  OF  THE  LATTICE.  INITIALIZE  THE  ZERO  STAGE  MODIFIED  NORM 
DO  10  T=0,N-1 
XO(T)  =  2.*S(T) 
TAU  =  TAU  +  S(T)**2 
10     CONTINUE 

DO  20  T=0,N 
IF(T.  EQ.  N)S(T)  =  0. 
IF(T.  EQ.  0)THEN 

X1(T)  =  S(T) 
ELSE 

X1(T)  =  S(T)  +  S(T-l) 
ENDIF 
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20     CONTINUE 
C    LOOP  TO  COMPUTE  THE  SYMMETRIC  PREDICTOR  COEFFICIENTS 
DO  101  K=1,M 
P(K,0)  =  1.0 
TTAU  =  TAU 
IF(K  . EQ.  KODD)THEN 

LL  =  (K+l)/2 
ELSE 

LL  =  K/2 
END  IF 
C     STORE  TAU(K-l),  AND  COMPUTE  TAU(K). 
TAU  =  0. 
DO  30  T=0,N+K-1 
TAU  =  TAU  +  X1(T)**2 
30     CONTINUE 

TAU  =  TAU/ 2. 
C     COMPUTE  ALPHA(K);  STORE  TAU(K)  FOR  NEXT  ITERATION. 
ALPHA(K)  =  TAU/TTAU 
TTAU  =  TAU 
C   LOOP  TO  CALCULATE  SINGULAR  PREDICTOR  COEFFICIENTS 
DO  40  1=1, LL 

P(K+1,I)  =  P(K,I)  +  P(K,I-1)  -  ALPHA(K)*P(K-1,I-1) 
C   DECISION  PATH  TO  CALCULATE  SYMMETRIC  SINGULAR  PREDICTOR  COEFFICIENTS 
IF(I.LT.  LL  .OR.  I.EQ.  K)GOTO  40 
L  =  K+l 
DO  50  J=I+1,K 
P(L,J)  =  P(L,L-J) 
50       CONTINUE 
40     CONTINUE 
C     LOOP  TO  COMPUTE  THE  (K+l)  PREDICTION  ERROR  VECTOR. 
DO  60  T=0,N+K 
IF(T  .EQ.  0)THEN 
AT(T)  =  X1(T) 
ELSEIF(T.  EQ.  N+K)THEN 
AT(T)  =  Xl(T-l) 
ELSE 

AT(T)  =  X1(T)  +  Xl(T-l)  -  ALPHA(K)*XO(T-l) 
ENDIF 
60     CONTINUE 
C     LOOPS  TO  UPDATE  PREDICTION  ERROR  VECTORS  FOR  NEXT  ITERATION. 
C         (SHIFT  XI  INTO  XO  AND  AT  INTO  XI) 
DO  70  T=0,N+K-1 
XO(T)  =  X1(T) 
70     CONTINUE 

DO  80  T=0,N+K 
X1(T)  =  AT(T) 
80     CONTINUE 

IF(K  .EQ.  KODD)THEN 

KODD  =  KODD  +  2 
ELSE 
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KEVEN  =  KEVEN  +  2 
ENDIF 
101    CONTINUE 
RETURN 
END 

SUBROUTINE  ASLAT( S ,M,N,PS , ALPHAS ,X1 ,AT,XO) 
C        SUBROUTINE  TO  COMPUTE  THE  LATTICE  REFLECTION  COEFFICIENTS 
C        USING  THE  ASYMMETRIC  PREDICTION  ERROR  VECTOR. 

REAL  X1(0:M+N),XO(0:M+N),PS(0:  100,0:  100)  ,  S(  0:  N)  ,  ALPHAS(M) 
REAL  AT(0:M+N) 
INTEGER  T 
C        INITIAL  CONDITIONS 
PS(0,0)  =  2. 
PS(1,1)  =  1. 
KODD  =  1 
KEVEN  =  2 
TAU  =  0. 
C   INITIALIZE  THE  PREDICTION  ERROR  VECTORS  FOR  THE  ZERO  AND  1ST 
C   STAGES  OF  THE  LATTICE.  INITIALIZE  THE  ZERO  STAGE  MODIFIED  NORM 
DO  10  T=0,N-1 
XO(T)  =  0. 

TAU  =  TAU  +  S(T)**2 
10     CONTINUE 

DO  20  T=0,N 
IF(T.  EQ.  N)X1(T)  =  -S(T-l) 
IF(T.  EQ.  0)THEN 

X1(T)  =  S(T) 
ELSE 

X1(T)  =  S(T)  -  S(T-l) 
ENDIF 
20     CONTINUE 
C    LOOP  TO  COMPUTE  THE  REFLECTION  COEFFICIENTS 
DO  101  K=1,M 
PS(K,0)  =  1. 
TTAU  =  TAU 
IF(K  .EQ.  KODD)THEN 

LL  =  (K-l)/2 
ELSE 

LL  =  K/2 
ENDIF 
C     STORE  TAU(K-l),  AND  COMPUTE  TAU(K). 
TAU  =  0. 
DO  30  T=0,N+K-1 
TAU  =  TAU  +  X1(T)**2 
30     CONTINUE 

TAU  =  TAU/ 2. 
C     COMPUTE  ALPHA(K)  ;  STORE  TAU(K)  FOR  NEXT  ITERATION. 
ALPHAS(K)  =  TAU/TTAU 
TTAU  =  TAU 
C   LOOP  TO  CALCULATE  SINGULAR  PREDICTOR  COEFFICIENTS 
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DO  40  1=1, LL 

PS(K+1,I)  =  PS(K,I)  +  PS(K.I-l)  -  ALPHAS(K)*PS(K-1,I-1) 
C   DECISION  PATH  TO  CALCULATE  SYMMETRIC  SINGULAR  PREDICTOR  COEFFICIENTS 
IF(I.LT.  LL  .OR.  I.EQ.  K)GOTO  40 
L  =  K+l 
DO  50  J=I+1,K 
PS(L,J)  =  -PS(L.L-J) 
50       CONTINUE 
40     CONTINUE 
C     LOOP  TO  COMPUTE  THE  (K+l)  PREDICTION  ERROR  VECTOR. 
DO  60  T=0,N+K 
IF(T  . EQ.  0)THEN 

AT(T)  =  X1(T) 
ELSEIF(T.  EQ.  N+K)THEN 
AT(T)  =  Xl(T-l) 
ELSE 

AT(T)  =  X1(T)  +  Xl(T-l)  -  ALPHAS(K)*X0(T-1) 
ENDIF 
60     CONTINUE 
C     LOOPS  TO  UPDATE  PREDICTION  ERROR  VECTORS  FOR  NEXT  ITERATION. 
C         (SHIFT  XI  INTO  XO  AND  AT  INTO  XI) 
DO  70  T=0,N+K-1 
XO(T)  =  X1(T) 
70     CONTINUE 

DO  80  T=0,N+K 
X1(T)  =  AT(T) 
80     CONTINUE 

IF(K  .EQ.  KODD)THEN 

KODD  =  KODD  +  2 
ELSE 

KEVEN  =  KEVEN  +  2 
ENDIF 
101     CONTINUE 
RETURN 
END 
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